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I .  INTRODUCTION 


Since  an  ELF  signal  from  a  remote  transmitter  is  received  over  a  wide  range  of 
azimuth  angles,  lateral  ionospheric  gradients  produced  by,  for  example, 
nuclear  depressions  can  produce  significant  effects  on  propagation  in  the 
lower  ELF  band.  This  is  because  the  Fresnel  zone  size  can  be  comparable  to 
the  transverse  dimensions  over  which  the  disturbed  ionosphere  changes  sig¬ 
nificantly.  It  has  been  common  practice  to  estimate  off -path  effects  by  using 
a  simple  surface -propagation  model  introduced  bv  Wait  [1964],  which  is  predi¬ 
cated  on  the  assumption  that  the  field  can  be  separated  into  lateral  and 
he  ight  -  dependent  factors.  A  number  of  workers  have  extended  and  used  this 
model  to  predict  the  vertical  electric  field  disturbance  produced  by  lateral 
gradients  [Field  and  Joiner  1979,  1982;  Pappert,  1980;  Ferguson,  Hitney  and 
Pappert,  1982;  Shellman,  1983],  These  studies  involve  either  the  numerical 
solution  of  an  integral  equation  or  numerical  solution  of  the  two-dimensional 
wave  equation  for  the  vertical  E  field.  The  methods  are  quite  general  in  the 
sense  that  rather  arbitrarily  shaped  disturbances  can  be  modeled.  This  study 
is  restricted  to  cy  1  indr ic a  1  ly  symmetric  disturbances,  and  the  formalism  of 
G reifinger  and  Greifinger  [1977]  is  explored  for  the  purpose  of  calculating 
both  the  vertical  electric  and  the  horizontal  magnetic  field  components  of  the 
RF  wave.  This  is  of  significance  in  ELF  studies  since  it  is  generally  the  RF 
magnetic  field  component  normal  to  the  great  circle  path  joining  transmitter 
and  receiver  which  is  measured.  The  Greifingers'  formulation  is  in  terms  of  a 
s ingle  -  component  vector  potential  from  which  all  field  components  are  obtained 
as  derivatives.  The  vector  potential  is  expressible  in  terms  of  an  infinite 
sum  of  products  of  cylindrical  and  harmonic  functions,  and  so  derivatives  are 
also  expressible  as  infinite  sums  of  well-known  functions.  A  major  purpose  of 
this  study  is  to  examine  the  practicality  of  calculating  those  sums  for  a 
variety  of  conditions. 
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Conmon  to  all  of  the  cited  works,  as  well  as  to  the  present  tudy ,  is  the 
neglect  of  excitation  factor  and  other  height-gain  effects  so  that  disturbance 
effects  reflect  solely  the  laterally  dependent  part  of  the  fields.  Common 
also  to  all  of  the  works  is  the  assumption  that  the  single  nonevane scent  ELF 
mode  is  pure  TM .  This  is  a  reasonable  approximation  under  ambient  and 
depressed  ionospheric  conditions  but  certainly  can  be  more  serious  for 
sporadic-E  type  environments  [Pappert  and  Moler,  1978],  The  current  study  is 
also  based  on  a  flat-earth  formulation  ana  therefore  is  restricted  to  ranges  < 

4  Mm,  where,  at  the  outer  limit,  errors  of  the  order  of  a  few  tenths  of  a  dB 
would  be  anticipated.  Shellman's  [1983]  work  is  unique  in  that  it  does  make 
full  allowance  for  spherical  spreading. 


Essential  theoretical  background  is  given  in  the  following  section.  In  par¬ 
ticular,  radial  variations  of  the  cy  1  indrically  symmetric  disturbances  are 
approximated  by  a  slab  model.  In  this  respect  the  present  program  is  a 
generalization  of  the  program  used  to  generate  the  results  in  Pappert  [1985]. 
A  shortcoming  of  the  program  is  that  C-raf'c  addition  theorem  for  Bessel  func¬ 
tions  [Watson,  1952],  which  forms  the  basis  of  the  calculation,  diverges  when 
the  distances  from  the  center  of  the  disturbance  to  transmitter  and  receiver 
are  equal.  For  some  configurations,  series  convergence  for  extended  ranges 
about  those  points  can  be  difficult  or  impossible  to  achieve.  Nonetheless, 
the  program  can,  with  care,  be  used  for  a  wide  variety  of  conditions  to  gain 
insights  into  the  effects  of  lateral  gradients  of  ELF  propagation. 


Program  input  and  output  is  described  in  sections  111  and  IV,  respectively. 
The  program  requires  eigenangle  inputs  for  both  the  ambient  and  disturbed 
regions  of  the  guide.  These  must  be  supplied  from  a  waveguide  program  such  as 
that  of  Morfitt  and  Shellman  [1976]  .  Provision  is  made  for  plotting 
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amplitudes  and  phases  of  the  ratio  of  disturbed  fields  to  undisturbed  fields 
as  a  function  of  range  for  fixed  location  of  the  disturbance  or  for  fixed 
ranges  as  a  function  of  location  of  the  center  of  the  disturbance  along  a 
straight  line.  As  illustrations  of  possible  applications  of  the  present 
i  program,  results  for  a  sporadic-E  type  environment,  for  a  nuclear  burst,  and 

I  for  a  weak  solar  proton  event  are  given  in  section  V.  A  program  listing  is 

|  given  in  Appendix  A  and  the  plotting  program,  apart  from  the  Hewlett  Packard 

International  Standard  Plotting  Package  (HPISPP)  components,  is  given  in 
Appendix  B. 
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II.  SUMMARY  OF  EQUATIONS 


Consiste..v  with  the  simplified  surface  propagation  model,  as  developed  by  the 
Greifingers,  it  is  assumed  that  the  single  nonevanescent  modal  field  com¬ 
ponents  which  propagate  within  the  ear  th  -  iono  sphe  re  waveguide  at  ELF 
frequencies  s„e  derivable  from  a  vector  potential  which  has  only  a  vertical 
component  which  can  be  factored  into  a  lateral  and  a  he  ight  -  dependent  com¬ 
ponent.  In  the  following,  only  the  behavior  of  the  lateral  component  will  be 
dealt  with  so  that  excitation  factor  and  height-gain  effects  are  ignored.  In 
source-free  regions,  the  lateral  component,  g,  of  the  vector  potential 
satisfies  the  reduced  two-dimensional  wave  equation, 

V  g(p,<J>)+k  S  (p,<f>)g(p,<j>)  =  0.  (1) 


where 
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1  f-  m  +  ^  s. 

p  dp  {  op)  p 2  g t2 


(2) 


The  free -space  wave  number  is  denoted  by  k  and  S  in  the  sine  of  the  eigenangle 
characterizing  the  single  nonevanescent  ELF  mode.  The  polar  coordinates  p  and 
^  are  shown  in  the  scattering  diagram  of  Figure  1.  Ear  th  -  curvature  effects 
are  ignored,  so  that  Figure  1  represents  a  plan  view  in  a  flat-earth  geometry. 
In  the  present  work,  S  is  taken  to  be  a  function  of  p  only  so  that  the  distur¬ 
bance  is  cyl  indr  i  cal  lv  .symmetric.  A  slab  approximation  will  be  vised  to  model 
the  radial,  or  p,  variation.  Thus,  the  disturbance  will  be  modeled  as: 
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2  2 
S  =  S, 


2  2 
S  =  S, 


P  <  Pi 


Pi  <  P  <  P; 


Pn-1  <  P  ^  pn 


2  2 

S  =  S..  ,  Pw  o  <  P  <  Pm  n 
M-l  M-2  -  M-l 


S  =  SM  Vl  <  ' 


where  is  che  eigenvalue  for  the  n  slab,  M  is  the  total  number  of  slabs, 


and  SM  is  the  ambient  eigenvalue. 


The  transmitter  is  taken  to  be  a  horizontal  electric  dipole  located  at  the 

origin  of  Figure  1  and  oriented  at  an  angle  S  relative  to  the  x-axis.  The 

receiver  is  taken  to  be  located  along  the  x-axis  at  a  distance  r  from  the 

transmitter.  The  disturbance  is  centered  a  distance  r  from  the  transmitter. 

o 

The  incident  wave  associated  with  the  dipole  source  is  derivable  from  the 
later?1  potential  function  (a  time  dependence  exp(iwt)  is  assumed  throughout 
this  work) . 


g.  =  cos ( -  6 ) H  J  ;(kSTr), 


(  2 ) 

where  Hj  is  the  Hankel  function  of  the  second  kind  of  order  1  and  the  sub¬ 
script  T  on  S  denotes  the  slab  in  which  the  transmitter  is  located.  In  this 
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connection,  the  subsequent  development  allows  for  the  transmitter  to  be 
beneath  the  disturbance  or  outside  its  projection  on  the  ground. 

Once  the  lateral  potential  function,  g,  is  determined,  the  laterally  dependent 
parts  of  the  ground  magnetic  field  components  and  are  given  by 


H 

P 


H<l> 


_dg 
dp  ' 


(5) 


Also,  to  good  approximation,  the  laterally  dependent  part  of  the  ground  verti¬ 
cal  electric  field  is  given  by 


Ez  =  'iwSR&-  (6) 

where  the  subscript  R  on  S  signifies  the  slab  in  which  the  receiver  is  lo¬ 
cated.  The  subsequent  development  also  allows  for  the  receiver  to  be  beneath 
the  disturbance  or  outside  its  projection  on  the  ground. 

The  total  lateral  potential  is  determined  by  adding  to  g^  the  scattered  com¬ 
ponent  which  satisfies  Equation  (1)  and  by  requiring  E ^  and  to  be 

continuous  at  each  slab  boundary.  It  is  remarked  that  the  scattered  component 
of  g  represents  the  total  component  in  all  slabs  except  the  slab  containing 
the  transmitter. 

The  scattered  component  of  g  in  slabs  1,  n  ^  1,  M,  and  M  is  written  as  (M  is 
the  total  number  of  slabs) : 

i)  n  =  1 ,  P  <  P\ 
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ii)  n  f  1 ,M,  pnl  <  P  <  Pr 


=  7  fc  ^H'1  cosmi  +  s  ^n')  sinmi I  J  (kS  p) 
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+  7  (c/^cosmA  +  s/;n^sinm<j>  (kS  p) 

L -  1  hm  T  hm  T  m  n 


iii)  n  =  M,  p  <  p 
m 


=  H_  (cim)cosm<i>  +  sin)sinm,*,)Hr)(kSn^)' 


The  c(n)'s  ar.d  s^n^'s  are  constants  to  be  determined  by  the  boundary  condi - 
jm  jm 

(  2  ) 

tions  and  the  J  's  (H  's)  are  Bessel  functions  (Hankel  functions  of  the 
m  m 

second  kind).  The  sum  over  m  runs  from  zero  to  infinity.  Continuity  of  the 
ratio  E^/H^  at  the  slab  boundaries  yields  (T  is  the  slab  containing  the  trans¬ 
mitter,  and  M  is  the  outermost  slab): 
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p  =  S  J  (ks  p  , ) fj ' (ks  ,p  ,)  +  7(n'1)H(2) ' (kS  ,P  n)l 
F  n  m  nFn-l'L  m  n-lFn-l  'm  m  n-rn-kj 


q  -  S  ,  J '  (kS  P  .)[J  (kS  ,p  .)  +  7(nl)H^2)(kS  lP  .  )1 
n  n  - 1  mv  n  n- 1  |_  m  n-lFn-l7  'm  m  v  n-lFn-l'J 


u  =  S  ,H:2)'(kS  p  ,)rj  (ks  .p  .)  +  7(n'1)H(2) (ks  ,p  . )1 
n-1  m  nFn-l  L  m  n-lFn-l'  'm  m  v  n-r  n-1  J 


v  -  S  H(2)(kS  p  ,  )  ("J '  (kS  ,p  .  )  +  7(n'1;)H(2)'(kS  ,»  ,  )1  . 

n  tn  n  n- 1  \_  m  n-lFn-l'  'm  m  v  n-lFn-l'J 


The  primes  occurring  in  the  preceding  equations  represent  derivatives  with 
respect  to  the  arguments. 


For  n  >  T,  continuity  o£  the  ratio  of  Ez/Hi  at  the  slab  boundaries  yields 


.  (M- 1 ) 


,  (M- 1 )  (M-l) 

'im  _  j_m 
(  M  - 1 )  “  .(M-l) 
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SM  ,H(2)'(kSMpM  1)H(2)(kSM  ,pM  .  )-SmH(2)  (kSMp,.  ^^'(kS.,  ,  pM  ,) 
M-l  m  M  M-l  m  M-l  M-l  Mm  MKM-1  m  M-l  M-l 

SmH(  )(kSMpM  , )J'(kSM  ,pM  .)  -  SM  ,H(  ^ ' (kS  p  1 )JM(kSM  .pM  .) 

M  m  M  M-l  m  M-l  M-l  M-l  m  MrM-l  M  M-l  M-l 


,  1<T<M-2(16) 


I 

% 


pi  -.qi 

ul  -  vl! 


l<T<M-3,  T<n<M- 
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where 


pi  =  S  H^CkS  p  )  rr(n+1)J'  (kS  ,p  )  +  H^;'(kS  ,p  )] 
r  n  m  n  n  m  m  n+rn  m  n+1  n  J 


‘.•I 


ql  -  S  -H*'-'  '  (kS  p  )rr(n+1)J  (kS  ,p  )  +  H^CkS  .p  )1 
^  n+1  m  n  n  |_  m  mv  n+l^n'  m  v  n+r  n'J 


ul  =  S  ,  J '  (kS  p  )rr(n+1)J  (kS  .p  )  +  H'-^CkS  ,.p  )1 
n+1  m  n*n [  m  mx  n+rn'  m  n+r  n  J 


| 

V. 


vl  =  S  J  (kS  p  ) fr(n+1)J' (kS  ,.p  )  +  H^‘;'(ks  .p  )] 
n  nr  n  n  |_  m  nr  n+lpn  m  n+rn  J 


The  choice  of  defining  7  and  T  by  their  respective  ratios  has  been  found  to 
yield  the  best  numerical  stability  when  using  an  L-U  decomposition  algorithm 
for  solving  the  linear  equations  imposed  by  the  continuity  of  E g  and  Hi  at  the 


$ 

'f 


slab  boundaries  p^,  ^  and  p^ .  The  conditions  imposed  on  Equations  (10),  (11), 

(16),  and  (17)  show  that  for  M  =  2,  neither  7’s  nor  T's  need  be  calculated. 

( 1 )  ( 2  ) 

In  that  case,  there  are  only  two  unknowns  (c:  ,  c  ,  )  associated  with  the 

J  jm  hm 

( 1 )  ( 2  ) 

cosine  terms  for  each  m  and  two  unknowns  (s:  ,  s7  ')  associated  with  the  sine 

j  m  hm 

terms  for  each  m.  Each  pair  in  this  instance  is  calculated  from  the  condition 
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The  continuity  equations  at  p ^  ^  and  are 


of  continuity  of  and  at  p j . 

established  by  noting  that  the  incident  lateral  potential  functions  for  radia¬ 
tion  at  an  angle  a  relative  to  the  x-axis  is: 


(2) 


gi  -  cos (a- S)Hj  (kSTr) 


(22) 


By  using  Graf's  theorem  [Watson,  1952],  g^  may  be  expanded  as  follows: 


i)  P  <  r 


g.  =  -2  cos  (e-6)/ _  e  '  (kS_r  )J  (kS_p)cosm<|> 

l  -  mm  1  o  m  1  1 


+  2  sin(e-S)/ _  mG  (kS  r  )J  (kS_p) sinmi 

-  m  I  o  m  1  1 

m 


/  rp  \  /  t  \ 

(p^  'cosm<j>  +  'sinmi) 
me  T  ‘ms  T 


(23) 


ii)  p  >  r 


=  -2  cos  ( e-6)/_  €mJm(kSTro)Hi  }  (kSTp)cosm<t> 
m 


+  2  sin (e-S)l mL^kS^r^H^  ^  (kS^,p) sinmtj> 


10 


\  /  'T  \  /  T  N 

[_  (qj, _  cosm<j>  +  q  sinm<|>), 


where 


0 . 5  m  =  o 


£m  1  m  /  o 


G'  (x>  =  Hv  '  (x)/x 
m  m 


L  (x)=J  (x)/x. 
mv  '  m  " 


It  is  remarked  that  Equations  (23)  and  (24)  diverge  for  rQ  =  P.  thereby 


restricting  the  generality  of  the  development  as  mentioned  in  the  introduc¬ 


tion.  By  using  these  decompositions  along  with  the  formula  for  given  in 


Equation  (5)  and  the  formula  for  Eg  given  in  Equation  (6),  the  continuity 


conditions  at  p T  ^  and  p T  yield  for  the  coefficients  of  the  cosine  terms: 


(T-l)u(2) 


(T)u(2) 


2sf.cos(e-5)e  H(2)'(kS„r  )J  (kSTp„  .) 
T  mm  T  o  m  T  T-l7 


ST - 1 C j m  [Jm(kST-lPT-l)+7m  ^m  (kST- l^T- 1  ^  +ST [Cjn/ Jm (kSTpT- 1 } +Chr/Hm  (kSTpT-  L}] 


.  (T) T , 


(T)u(2), 


2STcos(e-<5)€ X  ;'(kSTr  )J'(kS  p  ) 
T  mm  lorn  li-i 


V 

"1 


$ 

$ 
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ST[Cjm;Jm(kSTpT)+chm)Hm  }  (kSTpT)] 'ST+lchm  ^  fiS/  ^ Jm(kST+lpT)+Hm  ^  (kST+lpT)] 


2  (T+l) rr( 1+1) 


2S  cos (e -  5 )  £  J  '  (kS  r  )H^  )  (kS  p  ) 
I  m  m  T  o  m  TT 


ST [Cjm^ Jm(kSTpT^+Chn/Hm  (kbTpT}]  ' ST+lchm  ^  ;  Jm (kST+lpT) +Hm  ^  '  (kST+lpT)] 


,(T)U(2), 


(T+l) fp( T+l ) T , 


2S  cos (© - 5 ) €  J ' (kS  r  )Hk  1 ' (kSTp„) . 
f  m  m  T  o  m  T  T 


The  corresponding  equations  for  the  coefficients  of  the  sine  terms  are: 


(T-l)u(2) 


(T)u(2) 


-2St  sin(e-fi)mGm(kSTro)Jm(kSTpT_1) 


_ST-lsjm  } [Jm(kST-lpT-l)+7mT  1 }  Hm  } ' (kST- 1PT- 1}] +ST [S Jm(kSTpT- l)+Shm)Hm  )'(kSTPT-l)] 


,(T-1)U(2), 


=  -2S  sin(e-,5)mG  (kS  r  )J'(kS  p  ) 
1  m  Tom  TT-1 


ST[Sjm;Jm(kSTPT)+Shm)Hm  }  (kSTpT}] 'ST+lShm  U Jm(kST+lpT) +Hi  }  (kST+lpT}] 


2  (T+l) r_(T+l) 


•  2S_sin(e- 6 )mL  (kSTr  )H( * } (kS_pT) 
i  m  I  o  m  T  T 


ST[=J"7;+ST'>T+=^,Hr>' ■  +VT>] -ST+l=h»+1’  [C 


‘OOvOtf 


In  the  most  general  case  there  are  four  equations  and  four  unknowns.  Special 
cases  are : 


T 


1:  Use  Equations  (30)  and  (31)  for  cosine  coefficients  and  set 

0.  Use  Equations  (34)  and  (35)  for  sine  coefficients  and  set 

0 


( 1 ) 

T  -  2:  Set  7^  ’=  0. 

m 

T  =  M-l:  Set  r(M)=  0. 

m 

('M') 

T  ”  M:  Use  Equations  (28)  and  (29)  for  cosine  coefficients  and  set  c;  ™ 

jm 

0.  Use  Equations  (32)  and  (33)  for  sine  coefficients  and  set  s^M^=  0. 

It  should  be  observed  that  if  M  =  2  or  3,  two  of  the  special  cases  can  occur 
simultaneously.  For  example,  if  M  —  2  and  T  =  2  then  Equations  (28)  and  (29) 

( 1 ) 

would  be  used  for  the  cosine  coefficients  with  7  set  to  zero  and  Equations 

m  M 

( 1 ) 

(32)  and  (33)  would  be  used  for  the  sine  terms  with  7'  '  set  to  zero.  After 

m 

solving  for  the  unknown  coefficients  in  Equations  (28)  through  (31)  and  (32) 
through  (35),  we  obtain  all  remaining  unknown  coefficients  from  Equations 

(10),  (11),  (16),  and  (17).  With  the  cfn^'s,  c!' ' s ,  s^n^'s  and  s^n^'s  known, 

j  m  hm  j  m  hm 


the  scattered  part  of  the  lateral  potential  function  is  given  by  Equations  (7) 


through  (9).  The  latter  represent  the  total  potential  in  all  slabs  except 
that  containing  the  transmitter.  In  that  case,  the  total  potential  is  given 
by  the  sum  of  the  scattered  potential  and  the  primary  lateral  potential  func¬ 
tion  given  in  Equation  (4).  To  improve  convergence  properties  in  the 
neighborhood  of  rQ  =  p  for  slabs  n  y  nT>  it  is  preferable  to  write  Equations 

(7)  through  (9)  as: 

i)  n  =  1 ,  P  <  P i ,  n  y  n„ 


(c^}J  (kSlP) 
j  m  m  1 


(Ml 

me 

0) 

me 


)  cosmjzS 


( 1 ) 

+  ( s .  'j  (kS , p ) 
j  m  m  1 


( 1 ) 

ms 

O) 

ms 


)  sinm0 


+  cos(-5)H1(2)(kS1r) 


(36) 


iii)  n  =  M,  rM  <  p ,  n  ^  nT 


-L[ 

m  L 


(c^M)H(2)(kS  p)  - 

nm  m  M 


)  cosm^ 


+  (s^V^kS  p) 
hm  m  M 


^  )  sinmj*  +  cos(-7)Hj  )  (kSMr) 


where  the  p's  and  q's  are  defined  in  Equations  (23)  and  (24).  In  particular, 

the  p's  are  used  if  p  <  r  and  the  q's  if  p  >  r  . 
r  o  n  o 


The  RF  field  components  are  given  by  Equations  (5)  and  (6).  From  Figure  1  it 
will  ba  seen  that  the  horizontal  magnetic  field  components  and  are 

expressed  as  follows  in  terms  of  the  components  and  H^: 


Hr  =  H^cosy  +  H^siny 


H  =  -H.siny  +  H  cosy, 

y  0  P 


Amplitude  behavior  of  E  /E  U,  H  /H  U  and  H  /H  11  and  phase  behavior  of  their 

tr  2  2  y  y  IT 

reciprocals  are  the  principal  program  outputs,  where  the  superscript  u  sig¬ 
nifies  the  unperturbed  value. 


III.  DESCRIPTION  OF  INPUT 

Most  input  to  the  program  is  supplied  in  a  namelist  tile  of  16  or  fewer 
characters.  Examples  of  two  such  files  are  shown  in  Tables  1  and  2.  Table  1 
gives  input  for  a  range  variation,  and  Table  2  provides  an  example  for  a 
disturbance  variation.  The  namelist  inputs  are: 

rho(mxslab) . Cylindrical  slab  boundaries  in  km.  Also, 

mxslab  is  a  parameter  setting  the  dimension¬ 
ing  for  the  number  of  slabs.  The  number  or 
rho's  required  is  one  less  than  the  number 
of  slabs. 

iflag . Integer  flag.  Rui.6e  variation  is  given  by 

iflag  =  0,  disturbance  variation  by  iflag  f 

0. 

nrslab . Integer  value  of  number  of  cylindrical 

slabs . 

rmin . Starting  value  in  km  for  range  variation 

when  iflag  =  0.  When  iflag  ^  0  it  is  the 
range  value  in  km  for  which  calculations  are 
performed  as  the  disturbance  moves. 

rmax . Last  range  value  calculated  (in  km)  when 

iflag  =  0.  Not  used  otherwise. 


Table  1.  Sample  input  (range  variation) 


&datum 

rho-2750 . ,2800. ,2850. ,2900. ,2950. ,3000. ,3050. ,3100. , 

3150. . 3200. , 

3250. . 3300. .3350. .3400. .3450. .3500. .3550. , 

3600. . 3650. .3700. , 

3750. . 3800. , 
iflag-O, 
nrslab=23 , 
rmin-100 . , 
rmax-5000 . , 
rinc-50 . , 
smin-1802 . 78 . , 
smax=40 . , 
sinc-50 . , 
alpha-0 .  , 
deltal-0 .  , 
delta2-90 .  , 
cl-(l. ,0.), 
c2-(0. ,1.), 
xyint-3000 . , 

thta-(81. 305, -40.807) , (81 . 307 , -40 . 801) , (81 . 309 , -40 . 796) , (81.311,-40.784) , 
(81.317,-40.759) , (81.331, -40.708) , (81.358,-40.604) , 

(81.411,-40.395)  , 

(81.515,-39.992) , (81 . 703 , -  39 . 266)  , 

(82.007,-38.105) , (82.420,-36.560) , 

(82.861,-34.938) , (83.231,-33.611)  , 

(83.481,-32.727) ,(83.627,-32.218) , 

(83.704,-31.948) , ( 83 . 743 , - 31 . 812 ) , (83.763,-31.744) , (83 . 772 ,  -31 . 712) , 
(83.776,-31.696) , (83.778,-31.688) , (83 . 781 , -31 . 681)  , 
frqkhz-. 075 , 
skip-220. , 

&end 


mmm 
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Table  2.  Sample  input  (disturbance  variation) 


&datum 

rho=31 .25,62.5,93.75,125. , 156 . 25 , 187 . 5 , 218 . 75 , 250 . , 281 . 25 , 312 . 5 , 343 . 75 , 
375. ,406.25,437.5,468.75,500.  , 
iflag=l , 
nrslab=17 , 
rmin-1600 . , 
rmax=0 . , 
rinc=50 . , 
smin-- 1200 .  , 
smax=-600 . , 
sinc=100 . , 
alpha=- 57.44, 
deltal-0 . , 
delta2=90 . , 
cl-(l. ,0.) , 
c2-(0. ,1.), 
xyint=1250 .  , 

thta-( 59. 392805, -65. 552101) , (60.929764,-63.640642) ,(62.466723,-61.729183), 
(64.003682,-59.817724) , (65.540642,-57.906265) , (67.077601,-55.994806) , 
(68.61456, -54.083348) , (70.151519,-52.171889) , (71.688479, -50.26043) , 
(73.225438, -48.348971) , (74.762397,-46.437512) , (76.299356,-44.526054) , 
(77.836361,-42.614595) , (79.373275,-40.703136) , (80.910234,-38.791677) , 
(82.447193,-36.880218) , (83.984153, -34.96876) , 
frqkhz- . 075 , 
skip-50 . , 

&end 
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.Range  increment  in  km  when  if  lag  =  0.  Not 


used  otherwise. 


.When  iflag  /  0,  smin  is  the  starting  value 


in  km  of  disturbance  center  along  s  (see 


Figure  2).  When  iflag  =  0  it  is  the  value 


of  the  disturbance  center  along  s  for  which 


range  variations  are  performed. 


.When  iflag  =  0  smax  is  the  final  value  in  km 


of  disturbance  center  along  s  (see  Figure 


2).  Not  used  otherwise. 


Increment  in  km  by  which  disturbance  center 


is  moved  along  s  (see  Figure  2).  Not  used 


otherwise . 


alpha . Angle  in  degrees  denoting  slope  of  straight 


line  upon  which  the  disturbance  center  lies 


(see  Figure  2).  If  iflag  =  0,  the  center  is 


fixed.  If  iflag  ^  0,  the  center  moves  along 


s.  In  the  latter  case  the  direction  of 


motion  is  along  positive  s  if  alpha  >  0  and 


along  negative  s  if  alpha  <  0. 


deltal,  delta2 . The  program  allows  for  two  horizontal 


electric  dipole  sources  at  the  origin. 


Deltal  and  delta2  are  the  angles  in  degrees 


r-r>rio  hv  the  dipoles  with  the  x-axis  (see 


Figure  1).  Crossed  dipoles  could  have,  for 
example,  deltal  —  0°,  delta2  =  90°. 


cl,  c2 . Complex  amplitude  factors  of  the  two 

electric  dipole  sources.  Dipoles  90°  out  of 
phase  could  have,  for  example,  cl  =  (1. , 
0.),  c2  =  (0 . , 1 . ) . 

xyint . If  alpha  ^  0.,  xyint  is  the  intercept  of  the 

straight  line  s  with  the  x-axis.  In  Figure 
2  that  intercept  is  denoted  by  D.  If  alpha 
=  0°,  s  is  parallel  to  x,  and  xyint  is  the  y- 
intercept  of  the  path  s  .  A  given  value  of 
xO  and  yO  can  be  set  for  range  variation 
calculations  by  setting  alpha  =90°.  Then, 
xO  =  xyint  and  yO  =  smin. 

thta(mxslab) . Input  eigenangles  in  degrees  from,  for 

example,  a  waveguide  program  such  as  that  of 
Morfitt  and  Shellman  [1976],  An  eigenangle 
for  each  slab  is  requires: 


frqkhz . Radio  frequency  in  kHz. 

skip . Distance  in  km  about  r Q  =  P  (see  Figure  1), 


for  which  calculations  are  skipped  because 
of  poor  convergence  of  the  Bessel  and  Hankel 
function  expansions.  The  actual  distance 


skipped  along  s  if  iflag  f  0  or  along  r  if 


iflag  =  0  will  depend  upon  the  actual 
geometry . 


55? 


A  crucial  value  is  mmax ,  which  is  set  in  a  parameter  statement  in  the  sub¬ 
routines  mlinit,  mlcoef,  mlflds,  and  hj func .  It  controls  the  number  of  terms 
used  in  the  Bessel  and  Hankel  function  expansions.  In  examples  given  in  a 
later  section  of  this  report,  mmax  values  of  50  and  60  have  been  used. 
Generally,  increasing  mmax  allows  the  namelist  input  skip  to  be  reduced  only 
slightly,  because  convergence  of  the  Bessel  and  Hankel  function  expansion  is 
very  slow  in  the  neighborhood  of  p  =  r  .  For  some  variations  overflow 

problems  have  occurred  with  mmax  settings  =  80. 
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IV.  DESCRIPTION  OF  OUTPUT 


A  hard  copy  of  output  is  available  in  the  file  'mlout.'  It  is  particularly 
useful  when  concerned  about  sensitive  areas  of  calculation  which  may  criti¬ 
cally  depend  upon  the  input  skip  and  upon  the  parameter  setting  mmax.  Table  3 
shows  output  corresponding  to  the  input  of  Table  2.  The  first  section  of 
output  echoes  the  namelist  input  values.  Next  is  a  printout  of  mmax.  Each 
section  following  that  begins  with  the  slab  number  containing  the  transmitter 
(nt)  and  the  slab  number  containing  the  receiver  (nr) .  Following  that  are  the 
angles  y,  ®  and  0  in  radians  (see  Figure  1)  for  the  current  geometry.  Next 
comes  the  range  (r)  in  km,  the  distance  from  transmitter  to  center  of  distur¬ 
bance  (r  )  in  km  and  the  distance  from  the  center  of  the  disturbance  to  the 

receiver  (rhorcv  -  this  is  denoted  by  p  in  Figure  1)  in  km.  Next  come  values 
for  xO  and  yO  in  km,  which  are  the  x  and  y  coordinates  of  the  center  of  the 
disturbance  (see  Figure  1).  Next  the  location  of  the  center  of  the  distur¬ 
bance  along  the  s  axis  (scoord)  is  given  in  km.  Then  output  for  the 
disturbed- to  -  ambient  amplitude  ratios  in  dB  and  their  phase  differences  in 
degrees  are  given.  In  sensitive  regions  a  convergence  warning  is  given  and  a 
comparison  is  made  between  the  partial  sum  decomposition  of  the  fields  (ezps, 
hyps,  hrps)  and  their  exact  values  (eztst,  hytst,  hrtst)  .  These  comparisons 
are  made  only  as  an  indicator  of  possible  problems.  They  compare  the  partial 
sum  calculation  with  the  exact  value  for  the  range  r  using  the  transmitter 
eigenangle.  Generally,  inadequacies  of  the  calculations  for  particular 
ranges,  etc.,  are  indicated  by  erratic  behavior  of  the  plots  generated  by  the 
plotting  routine  discussed  subsequently.  Observe,  that  although  sine  =100., 
.scoord  jumps  from  -1000.  to  -700.  This  is  because  skip  =  50.  and  because  |r  - 

rhorcv|<50.  for  scoord  =  -900.  and  -800.  The  curves  generated  by  the  plotting 
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$ DATUM 

RHO  -  31.250000000000000,  62.500000000000000,  93.750000000000000, 

125.0,  156.250000000000000,  187.500000000000000,  218.750000000000000, 

250.0,  281.250000000000000,  312.500000000000000,  343.7500000000000 00, 

375.0,  406.250000000000000,  437.500000000000000,  468.750000000000000, 

500.0,  O.OOOOOOOOOOOOOOE+OOO,  O.OOOOOOOOOOOOOOE+OOO,  O.OOOOOOOOOOOOOOE+OOO, 
0 . 00000000000000E+000 , 

I FLAG  -  1, 

NRSLAB  -  17, 

RMIN  -  1600.0, 

RMAX  -  O.OOOOOOOOOOOOOOE+OOO, 

RINC  -  50.0, 

SMIN  -  -1200.0, 

SMAX  -  -600.0, 

SINC  -  100.0, 

ALPHA  -  -57.439999999999998, 

DELTAl  -  O.OOOOOOOOOOOOOOE+OOO, 

DELTA2  -  90.0, 

Cl-  (1.0, O.OOOOOOOOOOOOOOE+OOO) , 

C2  -  (O.OOOOOOOOOOOOOOE+OOO, 1.0) , 

XYINT  -  1250.0, 

THTA  -  (59.392805000000003,-65.552100999999993), 

(60 . 929763999999999 , -63 . 640642000000000) ,  (62 . 466723000000002 , 
-61.729182999999999) ,  (64.003681999999998, -59.817723999999998) , 
(65.540642000000005,-57.906264999999998) ,  (67.077601000000001, 
-55.994805999999997) ,  (68.614559999999997,-54.083348000000001) , 
(70.151518999999993,-52.171889000000000) ,  (71.688479000000001, 
-50.260429999999999) ,  (73.225437999999997,-48.348970999999999)  , 
(74.762396999999993,-46.437511999999998) ,  (76.299356000000003, 

-44 . 526054000000002) ,  (77 .836360999999997 , -42 . 614595000000001) , 

(79 . 373275000000007 , -40 . 703136000000001) ,  (80 . 910234000000003 , 
-38.791677000000000) ,  (82.447192999999999,-36.880217999999999) , 

(83 . 984153000000006 , -  34 . 968760000000003) ,  (0 . 00000000000000E+000 , 
O.OOOOOOOOOOOOOOE+OOO) ,  (O.OOOOOOOOOOOOOOE+OOO, O.OOOOOOOOOOOOOOE+OOO) , 

(0 . 00000000000000E+000 , 0 . OOOOOOOOOOOOOOE+OOO) , 

FRQKHZ  -  7 . 50000000000000E-002 , 

SKIP  -  50.0 
$end 

mmax”  5 1 
nt-  i7nr=  17 

chi-  . 793 15 75 7389 7 528 theta-  1 . 032309585380254phi-  1.316125581734791 

r-  1600. OrO-  1178 . 113912735650000rhorcv-  1419.356660749240700 

xO-  604. 180956552521020y0=  1011.393970280910480 

scoord-  -1200.0 

ez/ezu(db)-  -.165506486442835 

ez/ezu(deg)-  1.681470061692652 

hy/hyu(db)-  -3 . 94426540952422E-002 

hy/hyu(deg)=  2.882621844622661 

hr/hru(db)=  -.583301884765237 

hr/hru(deg)-  -3.085520568375234 

nt-  17nr-  17 

chi-  . 77 7432 2 56 2 65729 theta-  . 953569521683634phi-  1.410590963063211 

r-  1600. OrO-  1136 . 880831676240600rhorcv-  1321.703655468501400 

xO-  657 . 999210173144320y0=  927.111139424167960 

scoord-  -1100.0 

ez/ezu(db)-  -.252368993800051 

ez/ezu(deg)-  2.274333107647623 

hy/hyu(db)-  -6 . 70531978625971E-002 


hy/hyu(deg)-  3.439969724323339 
hr/hru(db)-  -.939223897961050 
hr/hru(deg)-  -3.454996754144538 
nt-  17nr-  17 

chi-  . 759203236734368theta=  . 869470069855074phi-  1.512919434423131 
r-  1600. OrO-  1103 . 197017528790900rhorcv-  1224.429571410443700 
xO-  711 . 817463793767620y0=  842.828308567425440 
convergence  warning 

eztst-  (.379935937456921, - ,415160527702075)ezps= 

( . 380245644175559 , - .416776157476905) 

hytst-  (-.270317510748229, . 384160851519301)hyps- 

(- .242660137752635, .370969437093575) 

hrtst-  (- .119815092087581, . 104044376392299)hrps= 

(- .106669605168189, .131705939060073) 
scoord-  -1000.0 
ez/ezu(db)=  -.349150339481386 
ez/ezu(deg)=  3.072486995179625 
hy/hyu(db)=  -.110864889350277 
hy/hyu(deg)-  4.108500334753247 
hr/hru(db)=  -1.342326254355762 
hr/hru(deg)=  -3.661623642820115 
nt-  17nr=  17 

chi=  . 68 19 132703 22 307 the ta=  . 594159466790211phi=  1.865520003900055 
r=  1600 . 0r0-  1053 . 888306054818000rhorcv=  936.060597793248460 
xO-  873 . 272224655637500y0-  589.979815997197870 
convergence  warning 

eztst-  (.379935937456921, - . 415160527702075)ezps= 

( . 379354113228671, - .415575175498021) 

hytst-  (- .270317510748229, . 384160851519301)hyps= 

(-.271096227011937, .401539362971401) 

hrtst-  (- .119815092087581, . 104044376392299)hrps- 

(- .137190518301539, .103207986744807) 

scoord-  -700.0 

ez/ezu(db)=  -.688823339745501 
ez/ezu(deg)-  6.319848973352530 
hy/hyu(db)-  -.405508957288416 
hy/hyu(deg)-  6.410720265774057 
hr/hru(db)=  -2.577415646978328 
hr/hru(deg)=  -4.538105647620277 
nt-  17nr=  17 

chi-  . 644465709277133 theta-  . 499356090273509phi=  1.997770941461931 
r-  1600. 0r0=  1056 . 042705429402500rhorcv-  841.746199995353440 
xO-  927 . 090478276260800y0=  505.696985140455350 
scoord-  -600.0 

ez/ezu(db)—  -.821916249843746 
ez/ezu(deg)-  7.486713866467464 
hy/hyu(db)-  -.588974755477543 
hy/hyu(deg)=  7.033520953575946 
hr/hru(db)-  -2.858686768921191 
hr/hru(deg)=  -5.432117372794706 


routine  will  be  plotted  by  using  simple  linear  interpolation  between  calcu¬ 
lated  points.  Thus  the  values  corresponding  to  scoord  =  -1000  and  scoord  =  - 
700  will  be  joined  by  a  straight  line. 


Another  output  of  the  program  is  an  unformatted  file  termed  "mlplot."  This  is 
used  in  conjunction  with  a  plotting  routine,  which  apart  from  components  of 
HPISPP  (Hewlett  Packard  International  Standard  Plotting  Package)  is  given  in 
Appendix  B  to  plot  any  combination  of  disturbed-over-ambient  amplitude  ratio 
or  phase  differences  for  either  range  or  source  variation.  The  lateral  field 
amplitude  ratios  are  in  dB ,  the  phase  differences  are  in  degrees,  and  the 
range  or  source  variation  is  in  Mm.  When  executing  the  plot  routine,  it  will 
first  request: 

Enter  name  of  file  containing  data  to  be  plotted. 

After  mlplot  is  entered,  the  screen  will  show: 

Do  you  want  amplitude  plots? 

Response :  y  or  n 

If  the  response  is  n,  the  questions  will  turn  to  phase  plots.  Otherwise 
questioning  for  amplitude  plots  will  continue  as  follows: 

How  many  amplitude  plots  do  you  want? 

Response:  1,  2  or  3 


1  ez  plot 

2  hy  plot 

3  hr  plot 


if  previous 
response 
was  1 


if  previous 
response 
was  2 


1  ez  and  hy  plots 

2  ez  and  hr  plots 

3  hy  and  hr  plots 


After  the  appropriate  code  (i.e.,  1,  2  or  3)  is  entered,  the  screen  will  show 

(size  x  and  size  y  are  dimensions  of  x-  and  y-axis  in  inches) : 

Current  values  for  size  x  and  size  y  are  6.0  6.0 
Do  you  want  to  change  them? 

Response :  y  or  n 

If  the  response  is  y  the  screen  will  read: 

Enter  values  for  size  x,  size  y 

After  size  x  and  size  y  are  entered,  or  if  response  to  previous  question  was 
n,  the  screen  will  then  read  (xmin  and  xmax  are  the  extreme  left  and  right  x 
coordinate  values ,  respectively) : 

Current  values  for  xmin  and  xmax  are  -5.0  7.0 
Do  you  want  to  change  them? 

Response:  y  or  n 

If  the  response  is  y  the  screen  will  read: 

Enter  values  for  xmin,  xmax 


a 

;y 


After  xmin  and  xmax  are  entered,  or  if  response  to  previous  question  was  n, 
the  screen  will  then  read  (ymin  and  ymax  are  the  extreme  bottom  and  top  y 
coordinate  values,  respectively): 


: 

: 

M 


Current  values  for  ymin  and  ymax  are  -10.0  4.0 
Do  you  want  to  change  them? 

Response  y  or  n 


If  the  response  is  y,  the  screen  will  read: 


$ 

v« 
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r 


.V* 


Enter  value  for  ymin,  ymax 


After  ymin  and  ymax  are  entered,  or  if  response  to  previous  question  was  n, 
the  screen  will  then  read: 


Tic  marks  will  be  every  2.0  units  on  the  x  axis  and  every  2.0  units  on 
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the  y-axis 

Do  you  want  to  change  them? 
Response  y  or  n 


I 


If  the  response  is  y,  the  screen  will  show: 


**; 

y, 

A 

5: 

y 


y 

y 

y 


Enter  values  for  xtic  and  ytic 


After  xtic  and  ytic  are  entered,  or  if  response  to  previous  question  was  n, 
amplitude  plots  will  then  be  started.  Next  inquiries  about  phase  plots  will 
appear  beginning  with 


5 

v 
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Do  you  want  phase  plots? 


L 


Response  y  or  n 

From  this  point  on,  the  interrogation  follows  the  same  sequence  described 
above  for  the  amplitude  plots.  It  asks  how  many  phase  plots  are  wanted,  which 
ones,  size  of  the  x-  and  y-axes,  coordinate  bounds,  and  tic  mark  locations. 
With  the  exception  of  the  slab  variation  sporadic-E  results,  the  plots  shown 
in  the  sample  results  section  of  this  report  were  produced  with  the  plotting 
package  described  above. 


V .  PROGRAM  LAYOUT 


This  section  describes  only  the  basic  features  of  the  ELF  propagation  program. 
In  particular,  the  plotting  routine  is  not  described  beyond  the  discussion  in 
the  previous  section. 

The  small  driver  program  main  opens  the  files  mlout  and  mplot  and  controls  the 
program  flow.  Subroutines  in  the  order  of  their  call  are  described  below: 

subroutine  mlinit 


Called  from  main,  mlinit  reads  in  namelist  data.  Quantities  independent  of 
transmitter  and  receiver  locations  are  calculated  in  mlinit.  These  include 

the  ' s  of  Equations  (10)  and  (11)  and  the  T^n^'s  of  Equations  (16)  and 

(17).  Note  that  the  r^n^'s  are  called  gammal(m,n)  in  the  program. 

subroutine  hifunc  (arg,  bj ,  bjp,  h2 ,  h2p) 

Bessel  functions,  bj  ,  and  Hankel  functions  of  the  second  kind,  h2 ,  and  their 
derivatives  bjp,  h2p  are  generated  in  hj  func  for  m  <  mmax  and  complex  argu¬ 
ment,  arg.  If  the  percentage  error  in  the  Wronskian  for  m  =  0  or  m  =  mmax 
exceeds  0.01%  a  warning  will  be  printed  in  mlout. 

subroutine  cbesiv  (z,  k,  bj ,  by,  kind,  nprint) 

Bessel  functions  of  the  first,  bj ,  and  second  kind,  by,  of  order  k  are  calcu¬ 
lated  by  standard  power  and  asymptotic  series  for  complex  argument,  z.  The 
routine  is  used  in  this  program  only  for  k  <  2  principally  for  the  purpose  of 
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obtaining  starting  values  and  as  checks  in  the  subroutine  hj  func .  h  the 
integer  flag  kind  =  1,  only  Bessel  functions  of  the  first  kind  will  be  gener¬ 
ated.  Any  other  value  of  kind  causes  generation  of  Bessel  functions  of  both 
kinds.  The  flag  nprint  =  0  causes  no  debug  printout,  whereas  nprint  —  1  will 
generate  debug  printout. 

subroutine  mlcoef 

In  this  subroutine  the  unknown  coefficients  in  Equations  (28)  through  (31)  and 
Equations  (32)  through  (35)  are  determined  by  a  Gaussian  elimination  process 

/  t  N  /Tii  \  /  T4.I  N 

implemented  in  clineq.  After  the  coefficients  c.  ,  c.  -  ,  c,  ,  c,  —  , 

r  j  m  j  m  hm  hm 

(T)  (T+l)  (T)  (T+l) 

s.  ,  s:  —  ,  s,  and  s,  -  have  been  determined  for  the  transmitter  slab  T 

j  m  j  m  hm  hm 

and  its  adjacent  slabs,  the  coefficients  for  slabs  n  f  T,  T+l  are  determined 
by  using  Equations  (10),  (11),  (16)  and  (17). 

subroutine  clineq  (a,  b,  x,  n,  ndim,  iflag,  err) 

The  solutions  of  simultaneous  linear  equations  with  complex  coefficients  are 
determined  by  clineq.  That  is,  it  solves  the  matrix  equation  a*x=t  for  the 
vector  x  of  length  n,  given  the  matrix  a  of  size  n  x  n  and  the  vector  b  of 
length  n  by  Crout’s  L-U  decomposition.  The  matrix  a  is  destroyed  by  clineq. 
The  integer  variable  ndim  must  always  be  greater  than  or  equal  to  n.  The 
integer  variable  iflag  has  been  set  to  zero  in  the  present  study.  The  real 
variable  err  is  computed  by  clineq  and  indicates  the  relative  error  in  the 
computed  solution  of  vector  x.  Printout  of  err  will  occur  only  if  it  exceeds 
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subroutine  mlflds 


All  quantities  necessary  for  calculating  the  field  components  are  available  in 
this  subroutine.  For  all  slabs  n  for  which  n  /  nT  the  total  lateral  potential 

g  equals  the  scattered  component  gs .  Thus  when  n  f  nT  the  field  components 
are  calculated  from  Equations  (5)  and  (6)  with  g  =  gg  and  g^  given  by 
Equations  (36)  through  (38).  When  n  =  n^,,  g  =  gg  +  g^,  where  g^  is  the  source 
field  given  by  Equation  (4).  Magnetic  field  components  and  are  calcu¬ 
lated  using  the  vector  resolution  given  by  Equations  (39)  and  (40).  Except 
for  namelist  and  mmax,  all  printout  onto  the  file  'mlout'  occurs  in  the  sub¬ 


routine  mlflds . 
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VI .  ILLUSTRATIVE  EXAMPLES 
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All  results  in  this  section  are  for  crossed  dipoles  90°  out  of  phase.  Figures 
3  and  4  repeat  results  with  the  present  program  of  results  at  75  Hz  obtained 
with  an  earlier  program  limited  to  two  slabs  [Pappert,  1985  ]  .  The  radius  of 
the  disturbance  is  500  km,  and  the  disturbance  models  a  resonant  sporadic-E 
environment  at  about  120  km  for  which  the  attenuation  rate  is  about  10  dB/Mm. 
The  complex  eigenangle  for  the  disturbed  region  is  (5Q.393°,  -65.553°)  and  the 
complex  eigenangle  for  the  ambient  region  is  (83.984°,  -34.969°).  Figure  3 
gives  amplitude  behavior  as  a  function  of  source  location  for  the  field  com¬ 
ponents  EZ ,  HY  and  HR  at  a  range  of  1.6  Mm  as  a  function  of  source,  location 
along  the  path  described  by  alpha  =  -57.44°  and  x  intercept  (D)  =  600  km  (see 
Figure  2).  Figure  4  shows  the  phase  behavior  for  the  same  disturbance  and 
path  of  travel.  The  deep  nulls  in  EZ  and  HY  and  their  large  phase  excursions 
have  been  reconciled  with  Wisconsin  test  facility  transmissions  to  Connecticut 
and  the  North  Atlantic  [Pappert,  1985].  The  amplitude  and  phase  behavior  of 
HR  does  not,  however,  appear  to  be  consistent  with  Bannister's  [1986]  radial 
magnetic  field  measurements.  The  reason  for  this  is  not  clear,  though  it  is 
suspected  that  it  may  be  the  result  of  TE  model  contamination  which  can  occur 
under  sporadic-E  environments  [Pappert  and  Moler,  1978], 


To  illustrate  convergence  with  slab  number  it  has  been  assumed  that  the 
sporadic-E  environment  of  the  previous  paragraph  is  such  that  the  eigenangles 
vary  linearly  with  radius  between  the  resonant  disturbed  value  at  the  origin 
and  the  ambient  value  at  the  radius  of  500  km.  Figures  5  through  10  show 
amplitude  and  phase  behaviors  of  all  field  components  for  slab  numbers  of  5, 
9,  and  17  as  a  function  of  source  locations  for  the  same  path  used  in  the 
previous  paragraph.  The  curves  indicate  reasonable  convergence  for  17  slabs. 
By  comparing  with  Figure  3,  it  will  be  seen  that  the  amplitude  fade  for  the 
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principal  field  components  EZ  and  HY  are  reduced  from  about  8  dB  to  between  2 
and  2.5  dB  by  virture  of  the  linear  variation  of  eigenangle.  Similarly, 
comparison  with  Figure  4  shows  that  the  phase  advance  for  EZ  and  HY  is  reduced 
from  more  than  60°  to  between  15°  and  20°  because  of  the  linear  variation. 

Further  examples  of  the  17-slab  model  of  the  linearly  varying  sporadic-E 


! 
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environments  of  the  previous  paragraph  are  shown  in  Figures  11  through  18. 


Figures  11  and  12  give  amplitude  and  phase  results,  respectively,  for  distur¬ 


bance  motion  along  a  path  which  passes  over  the  transmitter  so  that  the  slab 


index  containing  the  transmitter  changes.  There  is  some  evidence  of  this  in 


the  neighborhood  of  s  =  -500  km  and  s  =  100  km,  where  some  slight  discon¬ 


tinuities  are  visible.  Like  the  previous  results  the  EZ  and  HY  signatures  are 


remarkably  alike. 


Figures  13  and  14  show  amplitude  and  phase  results,  respectively,  for  distur¬ 


bance  motion  which  passes  over  the  receiver.  Unlike  previous  cases,  EZ  and  HY 


amplitude  and  phase  signatures  are  somewhat  dissimilar.  In  particular,  the  HY' 


amplitude  fade  is  significantly  deeper  than  the  EZ  fade,  and  the  EZ  fade 


advance  is  both  larger  and  broader  than  that  for  HY. 


Figures  15  and  16  give  amplitude  and  phase  results  as  a  function  of  range 


variation.  The  disturbance  center  is  located  at  xO  =  250  km  and  yO  =  0  km. 


It  is  believed  that  poor  convergence  of  the  partial  sum  expansions  is  respon¬ 


sible  for  the  kinks  which  appear  in  the  curves  in  the  neighborhood  of  700  km. 


In  this  connection  it  should  be  remarked  that  all  of  the  results  associated 


with  the  sporadic-E  environments  were  obtained  with  skip  =  50  km  and  mm  ax  - 
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Figures  17  and  18  also  show  amplitude  and  phase  behavior  as  a  function  of 
range.  The  disturbance  center  in  this  instance  is  located  at  xO  =  3750  km  and 
yO  =  0  km.  The  oscillatory  behavior  of  the  fades  in  the  region  between  the 
transmitter  and  the  disturbance  is  due  to  the  standing  wave  pattern  set  up  by 
the  incident  field  and  the  field  reflected  by  the  disturbance.  The  ambient 
waveguide  wavelength  is  3374  km,  and  it  will  be  seen  that  the  periodicity  of 
the  standing  wave  pattern  is,  as  expected,  about  one-half  of  that  value. 
Also,  the  EZ  and  HY  fades  are  dramatically  different  in  this  instance. 


To  illustrate  another  application  of  the  program,  Figures  19  through  26  show 
results  for  a  nuclear  environment  analyzed  by  the  Greifingers  (1977).  In 
particular  they  give  approximate  eigenvalues  at  45  Hz  for  a  bomb  fission  yield 
of  2  Mt  as  a  function  of  ground  truth,  assuming  that  the  fission  products  rose 
to  1000  km  at  10  min  after  burst.  The  slab  model  given  in  Table  4  is  based  on 
Table  3  given  in  the  Greifingers'  report.  In  this  connection,  eigenangles  are 
given  in  Table  4  rather  than  the  sine  of  eigenangles  given  by  the  Greifingers 
and  Table  4  has  been  expanded  ever  the  original  by  linearly  interpolating  the 
sine  of  the  eigenangle  between  inputs  in  Table  3  of  the  Greifingers'  report. 
The  original  Grief inger  data  would  allow  for  a  13- slab  model. 


Figures  19  and  20  shov  amplitude  and  phase  comparisons  of  the  principal  field 
components  EZ  and  HY  as  a  function  of  disturbance  motion  along  the  path 
described  by  o  =  -57.44°,  D  =  500  km  (see  Figure  2).  The  range  is  taken  to  be 
4  Mm .  In  these  calculations  mmax  =  60  and  skip  =  100  km.  The  maximum  EZ 
field  fade  is  about  3.5  dB ,  and  the  maximum  HY  field  fade  exceeds  4  dB ,  while 
phase  advances  approach  40°.  Jaggyness  due  to  slab  variation  is  evident.  The 
bulge  in  the  plots  close  to  S  =  3  Mm  is  indicative  of  partial  sum  convergence 
problems.  This  could  be  smoothed  out  by  choosing  a  larger  skip  value  (say  200 


mm. 
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Table  4.  Eigenangles  and  slab  radii  used  to  model  nuclear  environment 
sociated  with  2  Mt  of  fission  products  at  1000  km  altitude  10 
after  burst  time  (from  Greifinger  &  Greifinger,  1977). 


SLAB-N 

p (N) -Km 

Or(Deg) 

Qi (Deg) 

1 

400 

72.288 

-51.263 

2 

475 

72.356 

-51.085 

3 

550 

72.424 

-50.906 

4 

625 

72.493 

-50.726 

5 

700 

72.563 

-50.545 

6 

775 

72.759 

-50.308 

7 

850 

72.957 

-50.071 

8 

925 

73.157 

-49.832 

9 

1000 

73.359 

-49.593 

10 

1125 

73.192 

-49.204 

11 

1250 

73.020 

-48.813 

12 

1375 

72.845 

-48.419 

13 

1500 

72.667 

-48.024 

14 

1625 

73.274 

-47.591 

15 

1750 

73.895 

-47.160 

16 

1875 

74.529 

-46.731 

17 

2000 

75.178 

-46.305 

18 

2150 

75.351 

-45.906 

19 

2300 

75.528 

-45.502 

20 

2450 

75.707 

-45.095 

21 

2600 

75.889 

-44.683 

22 

2750 

76.283 

-44.361 

23 

2900 

76.683 

-44.039 

SLAB-N 

p (N) -Km 

9r (Deg) 

Qi(Deg) 

24 

3050 

77.090 

-43.718 

25 

3200 

77.504 

-43.396 

26 

3250 

77.978 

-42.701 

27 

3300 

78.468 

-41.995 

28 

3350 

78.976 

-41.279 

29 

3400 

79.501 

-40.553 

30 

3425 

79.746 

-40.098 

31 

3450 

79.996 

-39.637 

32 

3475 

80.251 

-39.172 

33 

3500 

80.512 

-38.700 

34 

3525 

80.832 

-38.411 

35 

3550 

81.156 

-38.120 

36 

3575 

81.487 

-37.830 

37 

3600 

81.823 

-37.540 

38 

3650 

81.842 

-37.121 

39 

3700 

81.861 

-36.696 

40 

3750 

81.880 

-36.266 

41 

3800 

81.897 

-35.831 

42 

3950 

82.057 

-34.910 

43 

4100 

82.222 

-33.961 

44 

4250 

82.391 

-32.981 

45 

4400 

82.565 

-31.969 

46 

00 

82.346 

-32.993 
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km).  It  could  in  principle  also  be  smoothed  out  by  choosing  a  larger  mmax , 
although  experience  has  shown  that  convergence  with  increasing  mmax  is  very 
slow . 

Figures  21  and  22  shows  EZ  and  HY  amplitude  and  phase  behavior  for  a  distur¬ 
bance  path  described  by  a  =  -57.44°,  D  =  3500  km  (see  Figure  2).  The  range  is 
4  Mm  so  that  the  center  of  the  disturbance  passes  within  about  420  km  of  the 
receiver.  Again  mmax  —  60  and  skip  =  100  km.  The  kink  in  the  HY  amplitude 
plot  and  the  break  in  the  derivatives  of  the  phase  plots  in  the  neighborhood 
of  S  =  -3  Mm  is  again  due  to  inadequacy  of  the  partial  sum  convergence. 
Jaggyness  due  to  slab  variation  is  also  evident. 

Figures  23  through  26  show  amplitude  and  phase  plots  as  a  function  of  range 
for  the  nuclear  depressed  environment.  For  Figures  23  and  24  the  disturbance 
is  centered  on  the  transmitter-receiver  path  at  a  distance  of  500  km  from  the 
transmitter.  For  Figures  25  and  26  the  disturbance  is  also  centered  on  the 
transmit  ter  -  rece  iver  path  but  at  a  distance  of  3500  km  from  the  transmitter. 
In  both  cases,  mmax  =  60  and  skip  =  100  km.  Amplitude  variations  fall  between 
about  1.5  and  -5  dB  and  phase  changes  between  about  -10°  and  40°. 

Considered  as  a  final  illustration  of  the  programs  application  is  a  polar  cap 
boundary  model  based  on  the  22  November  1982  weak  solar  proton  event  (SPE), 
whose  effects  on  ELF  propagation  were  examined  by  ray  trace  methods  by  Field, 
Warber,  and  Joiner  [1986],  Following  the  latter  reference,  the  diffuse  polar 
cap  boundary  is  assumed  to  extend  between  about  50°  and  55°  latitude,  and  to 
be  0 . 5  to  1  Mm  in  width,  and  the  sine  of  the  eigenangle  is  taken  to  be: 

S  =  S„pF  +  ( ^ amd  -  s„p_)/fl  +  expf -7 . 3 [ r-0 . 5(r i  +  r2)]/Ar)l,  (41) 


where 


S„nc,  =  1.25  -  0.117i  =  sine  of  eigenangle  in  disturbed  polar  cap  (42) 

Sr  E 

S.,,_  “  1.15  -  0.063i  =  sine  of  eigenangle  in  night  ambient  (43) 

Ari-D 

rj  =  3  Mm  (44) 

r2  =  rj  +  Ar  (45) 

Ar  =  transition  width  (46) 

and  r  is  measured  from  the  north  geomagnetic  pole. 

Figures  27  though  30  show  amplitude  and  phase  results  for  the  EZ  and  HY  com¬ 
ponents  for  a  hypothetical  motion  of  the  pole  along  a  path  described  by  a  = 
90°,  D  =  1802.78  km  (see  Figure  2).  The  calculations  have  been  performed  with 
mmax  «  60,  skip  =  50  km,  Ar  =  1  Mm,  and  33  slabs.  The  range  is  3  Mm  in 
Figures  27  and  28  and  4  Mm  in  Figures  29  and  30.  For  the  parameters  of 
Figures  27  through  30,  grazing  incidence  on  the  polar  cap  boundary  corresponds 
(the  case  of  interest  to  Field,  et  al .  [1986])  to  s  -  3  Mm.  It  will  be  seen 
that  the  HY  field  fade  there  is  about  1.2  dB.  This  corresponds  favorably  with 
the  ray  trace  estimate  of  1.4  dB  given  by  Field,  et  al .  [1986],  However,  the 
phase  advance  of  about  25°  to  30°  indicated  in  Figures  28  and  30  for  the 
grazing  incident  geometry  appears  to  be  inconsistent  with  measurements  which 
indicate  little  change  from  ambient. 

Figures  31  and  32  show  amplitude  and  phase  results,  respectively,  for  the  EZ 
and  HY  field  components  as  a  function  of  range.  Again  the  calculations  have 
been  performed  with  mmax  =  60,  skip  -  100  km,  Ar  -  1  Mm,  and  33  slabs.  The 
center  location  of  the  disturbance  is  given  by  xO  =  1802.78  km,  yO  —  3000  km. 
This  corresponds  to  grazing  incidence  on  the  polar  cap  boundary. 
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Irregularities  in  the  curves  in  the  neighborhood  of  the  range  =  3.6  Mm  is 
indicative  of  partial  sum  convergence  problems  ,  which  could  be  smoothed  out 
with  a  larger  value  of  the  input  parameter,  skip. 


Figures  33  through  36  are  analogues  of  Figures  27  through  30,  the  only  dif¬ 
ference  being  that  Ar  =  0.5  Mm  instead  of  1  Mm.  It  will  be  seen  that  the 
field  fades  for  5  =  3  Mm,  the  grazing  incidence  geometry,  is  <  1  dB ,  and  the 

phase  advance  is  somewhat  in  excess  of  25° .  This  is  to  be  contrasted  with  the 
ray  theory  results  of  Field,  et  al.,  which  suggest  6-dB  fades  for  the  polar 
cap  transition  thickness  Ar  =  500  km. 

Figures  37  and  38  are  the  analogues  of  Figures  31  and  32 .  They  have  been 
plotted  for  mmax  =  60,  skip  =  100  km,  number  of  slabs  —  33,  and  Ar  =  0.5  Mm. 
Again  the  evidence  would  point  to  fades  considerably  less  than  the  6  dB  sug¬ 
gested  by  ray  trace  arguments.  Again,  irregularities  in  the  neighborhood  of 
3.6  Mm  indicate  partial  sum  convergence  problems.  Other  wiggles  are  due  to 
slab-to-slab  variations. 


REFERENCES 


*  .'ii 


'•Vi 

'ft 

iSS 


*!• 

$: 

$ 

w 

»» 

s 


Ferguson,  J.A.,  L.R.  Hitney  and  R.A.  Pappert,  A  program  to  compute  vertical 
electric  ELF  fields  in  a  laterally  inhomogeneous  earth  ionosphere  waveguide, 
NOSC  TR  851,  prepared  for  the  Def.  Nucl .  Agency,  62  pp,  Dec  1982. 


Field,  E.C.,  C.R.  Warber  and  R.G.  Joiner,  "Focusing  and  shadowing  of  ELF 
signals,"  Radio  Sci . .  21  (3),  511-517,  1986. 


n 

M. 

A' 


Field,  E.C.  and  R.G.  Joiner,  "An  integral  equation  approach  to  long  wave 
propagation  in  a  non  stratified  earth  ionosphere  waveguide,"  Radio  Sci  ,  .  14 
(6),  1057-1068,  1979. 


Field,  E.C.  and  R.G.  Joiner,  "Effects  of  lateral  ionospheric  gradients  on  ELF 
propagation,"  Radio  Sci. .  17  (3),  693-700,  1982. 


m 

$ 


$ 

><!U 


'V  «l 

$ 

V) , 

& 

$ 

$ 


$ 

! 

i'.V 

><7J 

C 

liS 


Greifinger,  C.  and  P.  Greifinger,  Effects  of  a  cylindrically  symmetric  distur¬ 
bance  on  ELF  propagation  in  the  earth  ionosphere  waveguide,  Def.  Nucl.  Agency. 
Rep.  DNA  4399T,  40  pp.  R&D  Associates,  Marina  del  Rey,  CA,  June  1977. 

Morfitt,  D.G.  and  C.H.  Shellman,  'MODESRCH' ,  an  improved  computer  program  for 
obtaining  ELF/VLF/LF  mode  constants  in  an  earth  ionosphere  waveguide,  Interim 
report  77T,  prepared  for  the  Def.  Nucl.  Agency  by  the  Naval  Electronics  Lab. 
Cen.  (now  Naval  Ocean  System  Center),  227  pp,  Oct  1976. 


Pappert,  R.A.,  "Calculated  effects  of  traveling  spordic-E  on  nocturnal  ELF 
propagation.  Comparison  with  measurement,"  Radio  Sci .  .  20  (2),  229-246,  1985. 


Pappert,  R.A.,  "Effects  of  a  large  patch  of  sporadic-E  on  nighttime  propaga¬ 
tion  at  lower  elf,"  J.  Atmos.  Terr.  Phvs . .  42  (5),  417-425,  1980. 

Pappert,  R.A.  and  W.F.  Moler,  "A  theoretical  study  of  ELF  normal  mode  reflec¬ 
tion  and  absorption  produced  by  nighttime  ionospheres,"  J,  Atmos.  Terr.  Phvs . . 
40  (9),  1031-1045,  1978. 

She  lima1'.,  C.H.,  Propagation  of  ELF  waves  ir.  the  spherical  earth  -  ions  .-.phci'<- 
waveguide  with  transverse  and  longitudinal  variation  of  properties,  NOSC  TD 
607,  prepared  for  the  Def.  Nucl.  Agency,  194  pp .  June  1983. 

Wait,  J.R.  ,  "On  phase  changes  in  very  low-frequency  propagation  induced  by  an 
ionospheric  depression  of  finite  extent,"  J ,  Geonhvs ,  Res  .  .  59  (3),  441-445, 
1964. 

Watson,  G.H.  ,  A  treatise  on  the  theory  of  Bessel  functions.  Cambridge 


University  Press,  England,  1952,  p.  361, 


AMPLITUDE  RATIO-dB 


DELTA!®  0.00  DEG 
0ELTA2-  90.00  DEG 
X-INTERCEPT®  600.00  KM 


Cl»(  1.00  ,  0.00  ) 

C2®(  0.00  ,  1.00  ) 

ALPHA-  -57.44  DEG 

NRSLAB- 2 


4.0  |  T - 1 - 1 - 1 


DISTURBANCE  LOCATION-S  (Mm) 

I EZ/EZUI  _  FREQUENCY- 75.0  HZ 

I HY/HYUl  _  RANGE-  1.600  Mm 

I  HR/HRUI  . 


Figure  3.  Amplitude  ratios  vs.  disturbance  location  for  homogeneous  sporadic-E 
patch  of  500  km  radius. 
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Figure  4.  Phase  differences  vs,  disturbance  location  for  homogeneous  sporadic-E 
patch  of  500  km  radius. 
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Figure  5.  Slab  convergence  for  a  linearly  inhomogeneous  sporadic-E  patch  of  500 
km  radius,  EZ  amplitude  ratio. 
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Figure  6.  Slab  convergence  for  a  linearly  inhomogeneous  sporadic-E  patch  of  500 
km  radius,  EZ  phase  difference. 
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Figure  7.  Slab  convergence  for  a  linearly  inhomogeneous  sporadic-E  patch  of  500 
km  radius,  HY  amplitude  ratio. 
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Figure  8.  Slab  convergence  for  a  linearly  inhomogeneous  sporadic-E  patch  of  500 
km  radius,  HY  phase  difference. 
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Figure  9.  Slab  convergence  for  a  linearly  inhomogeneous  sporadic-E  patch  of  500 
km  radius,  HR  amplitude  ratio. 
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Figure  11.  Amplitude  ratios  vs.  disturbance  location  for  a  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  passes  over  the  transmitter. 
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Figure  12.  Phase  differences  vs.  disturbance  location  for  a  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  passes  over  the  transmitter. 
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Figure  14.  Phase  differences  vs.  disturbance  location  for  a  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  passes  over  the  receiver. 
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Figure  15.  Amplitude  ratios  vs.  r singe  for  an  on-path  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  overlaps  the  transmitter. 
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Figure  16.  Phase  differences  vs.  range  for  an  on-path  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  overlaps  the  transmitter. 
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Figure  18.  Phase  differences  vs.  range  for  an  on-path  linearly  inhomogeneous 
sporadic-E  patch  of  500  km  radius  which  overlaps  the  receiver. 
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Figure  21.  Amplitude  ratios  vs.  disturbance  location  for  a  nuclear  depression  of 
radius  4.4  Mm.  Disturbance  center  crosses  transmitter- receiver  path  at  3500  km. 
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Figure  22.  Phase  differences  vs.  disturbance  location  for  a  nuclear  depression  of 
radius  4.4  Mm.  Disturbance  center  crosses  transmitter-receiver  path  at  3500  km. 
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Figure  23.  Amplitude  ratios  vs.  range  for  a  nuclear  depression  of  radius  4.4  Mm 
Disturbance  center  is  on  path  and  500  km  from  transmitter. 
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Figure  27.  Amplitude  ratios  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  1  Mm,  range  =  3  Mm. 
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Figure  28.  Phase  differences  vs.  disturbance  location  for  a  weak  SPE.  Polar 
boundary  thickness  =  1  Mm,  range  =  3  Mm. 
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Figure  29.  Amplitude  ratios  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  1  Mm,  range  =  4  Mm. 
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Figure  30.  Phase  differences  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  1  Mm,  range  =  4  Mm. 
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Figure  31.  Amplitude  ratios  vs.  range  for  grazing  path  on  weak  SPE  polar  cap 
boundary  of  thickness  1  Mm. 
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Figure  32.  Phase  differences  vs.  range  for  grazing  path  on  weak  SPE  polar  cap 
boundary  of  thickness  1  Mm. 


WWW 


B 


Cl-(  1.00  .  0.00  ) 

C2- (  0.00  .  1.00  ) 

ALPHA-  90.00  DEG 
NRSLAB-  33 


0.0  r 


T 


DELTA1-  0.00  DEG 
DELTAS-  90.00  DEG 
X-INTERCEPT-  1802. 78KM 


~T 


"T" 


K 


CD 

"a 

i 

o 


< 

tr 


UJ 

a 


CL 

I 

< 


-4.0 


0.0 


DISTURBANCE  L0CATI0N-S  (Mm) 


4.0 


EZ/EZUI 

HY/HYUI 


FREQUENCY- 75.0 
RANGE-  3.000 


HZ 

Mm 


Figure  33.  Amplitude  ratios  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  0.5  Mm,  range  =  3  Mm. 
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Figure  34.  Phase  differences  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  0.5  Mm,  range  =  3  Mm. 
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Figure  35.  Amplitude  ratios  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  0.5  Mm,  range  =  4  Mm. 
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Figure  36.  Phase  differences  vs.  disturbance  location  for  a  weak  SPE.  Polar  cap 
boundary  thickness  =  0.5  Mm,  range  =  4  Mm. 
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Figure  37.  Amplitude  ratios  vs.  range  for  grazing  path  on  weak  SPE  polar  cap 
boundary  of  thickness  0.5  Mm. 
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Figure  38.  Phase  differences  vs.  range  for  grazing  path  of  weak  SPE  polar  cap 
boundary  of  thickness  0.5  Mm. 
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Figure  39.  Amplitude  ratios  vs.  range  for  penetrating  polar  cap  path. 
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APPENDIX  A 


ELF  SCATTERING  PROGRAM  LISTING 


main  routine  for  program 


implicit  real*8(a-h,o-z) 
open (unit-44,  file- 'mlout' ) 

open(unit— 40,file-'mlplot' , form- 'unformatted' ) 

call  mlinit 

call  mlcoef 

call  mlflds 

stop 


si 

M 


sfe 


yjy 

.!»© 

.vs* 

i 


subroutine  mlinit 
implicit  real*8(a-h,o-z) 

complex*16  cl , c2 , ampc , amps , thta , c , s , ssq ,ks , rllb , r22b , tl , 

$ lambda, £2 ,ksrho2 ,h2p ,h2 , bjp ,bj , j  rho2 , jprho2 ,hrho2 ,hprho2 , 

$ksrhol , j  rhol , jprhol , hrhol , hprhol ,num , den, gamma , dll , 

$ampprc , ampprs , gammal 
character*16  fn 
parameter(mmax-51) 
parameter(m''slab=20) 

dimension  siguia(mxslab)  ,  epsr  (mxslab)  ,rho(mxslab) 
dimension  c(mxslab) , s(mxslab) , thta(mxslab) , ssq(mxslab) , 

$ks (mxslab) , tl (mxslab) , lambda (mxslab) , f2 (mxslab) ,ksrho2 (mxslab) , 
$h2p(0 :ramax) ,h2(0:mmax) ,bjp(0:mmax) ,bj (0:mmax) ,jrho2(0:mmax, mxslab) 
$ ,jprho2(0 :mmax, mxslab) ,hrho2(0 :mmax, mxslab) ,hprhc2(0 :  inmax,  mxslab) , 
$ksrhol (mxslab) , j  rhol(0 :mmax .mxslab) , jprhol (0 :mmax .mxslab) , 

$hrhol(0 :mmax, mxslab) , hp rhol (0:mmax, mxslab) , gamma (0:mmax, mxslab) , 
$gammal(0:mmax, mxslab) 
common/one/ampc , amps ,ks , nrslab , rO , rho , s , ssq 

common/two/gamma .hprhol ,hprho2 , hrhol ,hrho2 .jprhol , jprho2 , j  rhol , 

$j  rho 2 , gammal 

common/ three/alpha , ampprc , ampprs , cl , c2 , cosalp , deltal , delta2 , 

$if lag, rinc , rmax, rmin, s inalp , sine , smax , smin.xyint , theta, xO ,y0 , skip 
data  pi/3. 1415926536/ 
data  vellit/2.997925e5/ 

namelist/datum/rho , if lag, nrslab , rmin, rmax , rinc , smin, smax , sine  , 
$alpha , deltal , delta2 , cl , c2 , xyint , thta , f rqkhz , skip 
do  n-l,nrslab 

read  in  slab  data 
end  do 
print  101 

format( ' input  filename:') 
read  102 , fn 
format(al6) 

open (unit-51 , f ile-fn, status-' old' ) 
read(51 , datum) 
write (44, da turn) 

write (40)  if lag, nrslab , f rqkhz , deltal , delta2 , cl ,  c2  ,  xyint ,  alpha 
f req-f rqkhz*1000 . 
waveno-2 . *pi*freq/vellit 
cosalp-cos(alpha*l . 7453292d-2) 
sinalp-sin(alpha*l .7453292d-2) 
if(alpha  .eq.  0.)then 
xO-smin 
yO-xyint 
else 

xO-smin*cosalp+xyint 
y0-smin*s inalp 
end  if 

rO-sqrt(xO**2+yO**2) 
if(r0  .It.  l.e-4)then 
if(x0  .It.  0.)then 
x0--l . e-4 

r0-sqrt(x0**2+y0**2) 

else 

xO-1 .e-4 

rO-sqrt (xO**2+yO**2 ) 
end  if 


r ' .li1  .»•* ,iil.<<Vi 


i 


‘5 


.« 

s 

s 


go  to  2 
end  if 

if(xO  .eq.  0.)then 
theta-pi/2 . 

if(yO  .It,  0 .) theta-- theta 
else 

theta-atan2 (yO ,x0) 
end  if 

write (44,*) ’xO-' ,x0 
write (44,*) 'yO-' ,y0 
deltal-deltal*l . 7453292e-2 
delta2-delta2*l . 7453292e-2 
ctmdl-cos (theta-deltal) 
c tmd2-cos ( theta- delta2) 
stmdl-s in (theta-deltal) 
s  tmd2-s in ( the  ta-delta2) 
s  the  ta-s in ( the  ta ) 
c  the  t a-c  o  s ( theta) 
betal--deltal 
beta2--delta2 
sbetal-sin(betal) 
sbeta2-sin(beta2) 
cbetal-cos(betal) 
cbeta2-cos(beta2) 
ampc-2 . *ctmdl*cl+2 . *ctmd2*c2 
amps-2 . *stmdl*cl+2 . *stmd2*c2 
ampprc-cl*cbetal+c2*cbeta2 
ampprs=-cl*sbetal-c2*sbeta2 
do  n-l,nrslab 

thta(n)— thta(n)*l . 7453292e-2 
c(n)-zcos(thta(n)) 
s(n)-zsin(thta(n) ) 
ssq(n)— s(n)*s(n) 
ks(n)-waveno*s (n) 

call  rbars(a(n) ,c(n) ,sigma(n) ,epsr(n) , freq , rllb , r22b) 
dll-(l ,+rllb)*(l ,+rllb) 
lambda (n) -- dll* tl (n)*s(n) 
f2(n)-c(n)*(l.-rllb)/(l.+rllb) 


if(n  .eq.  l)then 

ksrho2 (n)=ks (n)*rho(n) 

call  hjfunc(ksrho2(n) ,bj ,bjp,h2 ,h2p) 

do  m-0 , mmax - 1 

jrho2(m,n)=bj (m) 
jprho2(m,n)-bjp(m) 
hrho2(m,n)-h2(m) 
hprho2(m,n)-h2p(m) 
end  do 
else 

if(n  .eq.  nrslab)then 
ksrhol(n)-ks(n)*rho(n-l) 
call  hj func(ksrhol(n) ,bj , bjp ,h2 ,h2p) 
do  m-0 , mmax - 1 

jrhol(m,n)-bj (m) 
jprhol(m, n)-bjp(m) 
hrhol(m,n)-h2(m) 
hprhol(m,n)-h2p(m) 
end  do 


W3HWOT 


else 

ksrho2(n)-ks(n)*rho(n) 
ksrhol(n)-ks(n)*rho(n-l) 
call  hjfunc(ksrho2(n) ,bj ,bjp,h2 ,h2p) 
do  m-0 , mmax - 1 

jrho2(m,n)-bj (m) 
jprho2(m,n)-bjp(m) 
hrho2(m,n)-h2(m) 
hprho2(m,n)-h2p(m) 
end  do 

call  hj func(ksrhol (n) ,bj ,bjp,h2 ,h2p) 
do  m-0, mmax -1 

jrhol(m,n)-bj (m) 
jprhol(m,n)-bjp(m) 
hrhol (m , n)-h2 (m) 
hprhol  (m ,  n)-h2p  (in) 
end  do 
end  if 
end  if 


$ 

$ 

$ 

$ 


$ 


if(n  .ne.  1  .and.  n  ,ne.  nrslab)then 
if(n  .eq.  2)then 
do  m-0, mmax -1 

num-s(2)*jprho2(m, l)*jrhol (m, 2) 
num— num-s(l)*j rho2(m, l)*jprhol(m, 2) 
den-s ( 1 ) * j  rho2 (m , 1 ) *hprhol ( m , 2 ) 
den-den- s (2) *jprho2(m, l)*hrhol(m, 2) 
gamma ( m , n ) -num /den 
end  do 
else 

do  m-0, mmax -1 

num-s(n)*jrhol (m, n)*(jprho2 (m, n-l)+gamma(m,n-l) 
*hprho2(m,n-l)) 

num-num- s(n-l)*jprhol(m,n)*(jrho2(m,n-l)+gamma(m,n-l) 
*hrho2(m,n-l) ) 

den-s (n-l)*hprhol(m,n)*(jrho2(m,n-l)+gamma(m,n-l) 
*hrho2(m,n-l)) 

den-den- s (n) *hrhol (m , n) * ( j  prho2 (m , n- 1 ) +gamma (m , n - 1 ) * 
hprho2(m,n-l)) 
gamma ( m , n ) -num/ den 
end  do 
end  if 
end  if 
end  do 

do  n-nrslab-1 , 2 , -1 

if(n  .eq.  nrslab-l)then 
do  m-0 , mmax - 1 

den-s (nrs lab ) *hrho 1 (m, nr slab) *jprho2 (m, nrslab- 1) 
den-den- s (nrslab - 1 ) *hprhol (m , nr slab) *j rho2 (m , nrslab - 1 ) 
num-s (nrs lab  - 1 ) *hprhol (m , nrs lab ) *hrho2 (m , nr slab  - 1 ) 
num-num- s ( nrslab ) *hrhol (m , nr slab ) *hprho2 (m , nr slab  - 1 ) 
gammal(m,nrslab-l)-num/den 
end  do 
else 

do  m-0, mmax -1 

rum-s (n+1 ) * j  prho2 (m , n) * ( gammal (m , n+1 ) * j rhol (m , n+1 )  + 
hrhol (m,n+l) ) 

num-num- s(n)*j rho2(m,n)*(gammal(m,n+l)*jprhol (m,n+l)+ 
hprhol (m, n+1) ) 


$ 


den-s  (n)  *hrho2  (m , n) * ( gammal  (m ,  rt+1 )  *  j  prhol  (m ,  n+1 )  + 
hprhol(m,n+l) ) 

den-den- s (n+1 ) *hp rho2 (ra,n)*( gammal (m, n+1 )*j rhol (m,n+l)+ 
hrhol(m,n+l) ) 
gammal (m , n)-den/num 
end  do 
end  if 
end  do 
return 


subroutine  mlcoef 
implicit  complex*16 (a-h , o- z) 

complex*16  ks , j rO , j prO , IrO , j rhol , j  rho2 , j prhol , j prho2 , num 
real*8  rO,rho 
parameter (mmax=51) 
parameter (mxslab-20) 

dimension  rho(mxslab) ,ks(mxslab) , jr0(0:mmax) , jpr0(0:mmax) , 

$hr0(0 :mmax) ,hpr0(0 : mmax) , lr0(0 : mmax) , gr0(0 :mmax) , 
$jrhol(0:mmax,mxslab) , j prhol (0 : mmax, mxslab) , jrho2(0:mmax,mxslab) , 

$ j prho2 (0: mmax, mxs lab) ,hrhol(0: mmax, mxslab) , hprhol (0 : mmax , mxslab) , 
$hrho2(0: mmax, mxslab) ,hprho2(0 : mmax, mxs lab) , ssq (mxslab) , s (mxslab) , 
$a2 (2 , 2) ,b2 (2) ,x2 (2) , a4(4 , 4) ,b4 (4) , x4(4) , xj  c (0 : mmax , mxs lab) , 
$xjs(0: mmax, mxslab) ,xhc(0: mmax, mxslab) ,xhs (0 : mmax .mxslab) , 

$ gamma (0 ; mmax .mxslab) , a2sav(2 , 2) , a4sav(4 , 4) , gammal (0 : mmax , mxslab) 
common/one/ampc , amps , ks , nr slab , rO , rho , s , ssq 

common/ two/gamma , hprhol ,hprho2 ,hrhol ,hrho2 , j prhol , jprho2 , j  rhol , 

$j rho2 , gammal 

common/four/xhc ,xhs , xj  c , xj  s 
common/f ive/grO ,hpr0 ,hr0 , jprO , j  rO , IrO 
nt-nrslab 
do  n-1 ,nrslab-l 

if(r0  . le.  rho(n))then 
nt=n 
go  to  1 
end  if 
end  do 
continue 
arg-ks (nt)*r0 

call  hj  func(arg, jrO , jprO ,hr0 ,hpr0) 
do  m-0 , mmax - 1 

Ir0(m)-jr0(m) /arg 
gr0(m)-hr0(m)/arg 
end  do 

if(nt  .eq.  1  .and.  nrslab  .eq.  2)then 
do  m-0 , mmax - 1 

a2(l , l)-ssq(l)*j rho2(m, 1) 
a2<l , 2)=-ssq(2)*hrhol(m, 2) 
a2(2 , l)»s(l)*jprho2(m, 1) 
a2(2,2)--s(2)*hprhol(m,2) 
a2sav(l ,  l)>=a2(l ,  1) 
a2sav(l,2)-a2(l,2) 
a2sav(2 , l)-a2 (2 , 1) 
a2sav(2 , 2)-a2(2 ,2) 

b2(l)-ampc*ssq(l)*jprO(m)*hrho2(m, 1) 
b2(2)-ampc*s(l)*jpr0(m)*hprho2(m, 1) 
if(m  .eq.  0)then 
b2(l)-.5*b2(l) 
b2(2)-. 5*b2(2) 
end  if 

call  clineq(a2 ,b2 ,x2 , 2 , 2 , 0 , err) 

xjc(m,l)-x2(l) 

xhc(m, 2)-x2(2) 

b2(l)--amps*ssq(l)*m*lrO(m)*hrho2(m, 1) 
b2(2)--amps*s(l)*m*lrO(m)*hprho2(m, 1) 
call  clineq (a2sav,b2 , x2 , 2 , 2 , 0 , err) 
xjs(m, l)-x2(l) 
xhs(m, 2)-x2(2) 


end  do 
go  to  2 
end  if 

if(nt  .eq.  2  .and.  nrslab  .eq.  2)then 
do  m-0 , mmax - 1 

a2(l,l)=-ssq(l)*i rho2 (m, 1) 
a2 (1 , 2)=ssq(2)*hrhol(m, 2) 
a2 (2 , l)=-s(l)*jprho2(m, 1) 
a2(2,2)«s(2)*hprhol(m,2) 
a2sav(l , l)-a2 (1 , 1) 
a2sav(l , 2)=a2 (1 , 2) 
a2sav(2 , l)=a2 (2 ,1) 
a2sav(2,2)-a2(2,2) 

b2(l)-ssq(2)*ampc*hpr0(m)*j rhol(m, 2) 
b2(2)-s(2) *ampc*hprO (m) *j prhol (m , 2 ) 
if  (in  .eq.  0)then 
b2(l)=.5*b2(l) 
b2(2)=. 5*b2(2) 
end  if 

call  clineq(a2 , b2 , x2 , 2 , 2 , 0 , err) 
xjc(m, l)-x2(l) 
xhc (m, 2)“x2 (2) 

b2 (l)~-amps*ssq(2)*m*gr0(m)*j rhol (m, 2) 
b2(2)  =  -amps*s(2)*m*gr0(ni)*jprhol(m) 2) 
call  clineq(a2sav,b2 ,x2,2,2,0,err) 
xjs(m, l)-x2(l) 
xhs(m, 2)=x2(2) 
end  do 
go  to  2 
end  if 


if(nt  .eq.  2  .and.  nrslab  .eq.  3)then 
do  ra-O.mmax-l 

a4(l,l)=-ssq(l)*j rho2(m , 1) 
a4(l , 2)=ssq(2)*j rhol(m, 2) 
a4(l , 3)=ssq(2)*hrhol(m, 2) 
a4(l,4)-0. 

a4(2 , l)=-s(l)*jprho2 (m, 1) 
a4(2 , 2)-s (2 )*j prhol (m, 2) 
a4(2 , 3)-s(2)*hprhol(m, 2) 
a4(2,4)-0. 
a4(3,l)-0. 

a4(3 , 2)~ssq(2)*j rho2(m,2) 
a4( 3 , 3)«=ssq(2)*hrho2  (m,  2) 
a4(3 ,4)-- ssq(3)*hrhol(m, 3) 
a4(4,l)-0. 

a4(4 , 2)-s(2)*jprho2(m, 2) 
a4(4 , 3)-s (2)*hprho2 (m, 2) 
a4 ( 4 , 4 ) “ - s ( 3 ) *hprho 1 ( ra , 3) 
a4sav( 1 , l)-a4(l , 1) 
a4sav(l,2)-a4(l,2) 
a4sav( 1 , 3)-a4(l , 3) 
a4sav(l , 4)-a4(l ,4) 
a4sav(2 , l)«a4(2 , 1) 
a4sav(2 , 2)-a4(2 , 2) 
a4sav(2,3)-a4(2,3) 
a4sav(2 , 4)-a4(2 ,4) 
a4sav(3 , l)-a4(3 , 1) 
a4sav(3 , 2)-a4(3 , 2) 
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a4sav(3 , 3)-a4(3 , 3) 
a4sav(3 ,4)-a4(3 ,4) 
a4sav(4, l)-a4(4, 1) 
a4sav(4 , 2)-a4(4 ,2) 
a4sav(4 , 3)=-a4(4 , 3) 
a4sav(4 ,4)-a4(4,4) 

b4(l)-ampc*ssq(2)*hprO(m)*j rhol(m, 2) 
b4(2)-ampc*s(2)*hpr0(m)*jprhol (m, 2) 
b4 ( 3 )~ampc*ssq ( 2 ) * j  prO (m) *hrho2 (m , 2 ) 
b4(4)-ampc*s(2)*jpr0(m)*hprho2 (m , 2) 
if(m  .eq.  0)then 
b4(l)-. 5*b4(l) 
b4(2)-.5*b4(2) 
b4(3)-. 5*b4(3) 
b4(4)-.5*b4(4) 
end  if 

call  clineq(a4 ,b4, x4 , 4 , 4 , 0 , err) 
xjc(m, l)“x4(l) 
xjc(m,2)-x4(2) 
xhc(m, 2)-x4(3) 
xhc(m, 3)=*x4(4) 

b4(l)--amps*ssq(2)*m*gr0(m)*j rhol(m, 2) 
b4(2)«-amps*s (2)*m*grO(m)*jprhol(m, 2) 
b4(3)--amps*ssq(2)*m*lrO(n))*hrho2  (m,  2) 
b4(4)--amps*s (2)*m*lrO(m)*hprho2 (m, 2) 
call  clineq(a4sav,b4 ,x4 ,4 ,4 , 0 , err) 
xjs(ra,l)“x4(l) 
xjs(m,2)-x4(2) 
xhs(m,2)”x4(3) 
xhs(m, 3)-x4(4) 
end  do 
go  to  2 
end  if 

if(nt  .eq.  1  .and.  nrslab  .ge.  3) then 
do  m=0,mmax-l 

a2(l,l)=ssq(l)*j  rho2 (m, 1) 

a2 (1 , 2)“-ssq(2)*(gammal(m , 2)*j rhol(m , 2)+hrhol (m, 2) ) 
a2(2 , l)-s(l)*jprho2(m, 1) 

a2(2 , 2)--s(2)*(gammal(m, 2)*jprhol (m , 2)+hprhol(m, 2) ) 

a2sav(l , l)-a2(l , 1) 

a2sav(l , 2)-a2(l , 2) 

a2sav(2,l)-a2(2,l) 

a2sav(2,2)-a2(2,2) 

b2(l)-ampc*ssq(l)*jprO(m)*hrho2  (m,  1) 
b2(2)-ampc*s(l)*jpr0(m)*hprho2(m, 1) 
if(m  .eq.  0)then 
b2(l)-. 5*b2(l) 
b2(2)-.5*b2(2) 
end  if 

call  clineq(a2 , b2 , x2 , 2 , 2 , 0 , err) 

xjc(m,l)-x2(l) 

xhc(m, 2)-x2(2) 

b2(l)--amps*ssq(l)*in*lrO(ra)*hrho2(m,  1) 
b2(2)--amps*s(l)*m*lr0(ni)*hprho2  (m,  1) 
call  clineq(a2sav, b2 , x2 , 2 , 2 , 0 , err) 
xjs(m, l)-x2(l) 
xhs(m,2)-x2(2) 

xjc(m,  2)-gaumial(m,  2)*xhc(m,  2) 
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xj  s  (m,  2)-gammal(m,  2)*xhs  (m,  2) 
do  n=3 , nrslab- 1 

num«ssq(n-l)*(gammal(m ,n-l)*j rho2 (m , n- l)+hrho2 (m, n- 1) ) 
den-ssq(n)*(gammal(m, n)*j rhol(m,n)+hrhol(m,n) ) 
xhc  ( m ,  n )  =num*xhc  (m,n-l)/der> 
xjc(m,n)=gammal(m,n)*xhc(m,n) 
xhs (m , n)-num*xhs (m, n- l)/den 
xj  s  (m ,  r.)  -gammal  (m ,  n)  *xhs  (m ,  n) 
end  do 

num-ssq(nrslab - l)*(gammal (m, nrs] ab- l)*j  rho2 (m,nrslab-l)  + 
hrho2 (m, nr slab -1)) 
den-ssq(nrslab) *hrhol (m , nr slab) 
xhc (m , nr s lab ) =num*xhc (m,nrslab-l) /den 
xhs (m , nr s lab ) -num*xhs (m,nrslab-l) /den 
end  do 
go  to  2 
end  if 


if(nt  .eq.  nrslab  .and.  nrslab  . ge .  3) then 
do  m-0,mmax-l 

a2 (1 , 1)=- ssq (nrslab- 1)*( jrho2 (m, nr slab- l)+gamma(m, nrslab- 1)* 
hrho2(ra, nrslab- 1) ) 
a2(l , 2)-ssq(nrslab)*hrhol(m,nrslab) 

a2(2 , 1)=- s (nrslab -1)*( jprho2(m, nrslab -l)+gamma(m, nrslab -1)* 
hprho2(m,nrslab-l)) 
a2(2 , 2)=s(nrslab)*hprhol(m,nrslab) 
a2sav(l , l)-a2 (1 , 1) 
a2sav(l , 2)-a2(l , 2) 
a2sav(2 , l)=a2(2 , 1) 
a2sav(2,2)-a2(2,2) 

b2(l)-ssq(nrslab)*ampc*hpr0(m)*j rhol(m, nrslab) 
b2 (2 )-s (nrslab) *ampc*hprO(m)*jprhol(m, nrslab) 
if(m  .eq.  0)then 
b2(l)-.5*b2(l) 
b2(2)“. 5*b2(2) 
end  if 

call  clineq(a2 , b2 ,x2 , 2 , 2 ,0 , err) 
xjc(ra, nrslab- l)=x2(l) 
xhc (m,nrslab)-x2(2) 

b2(l)--ssq(nrslab)*amps*m*gr0(m)*j  rhol(m,nrslab) 
b2 (2 )~-s (nrslab) *amps*m*grO(m)*jprhol(m, nrslab) 
call  clineq(a2sav,b2 , x2 , 2 , 2 , 0 , err) 
xjs(m, nrslab- l)=x2(l) 
xhs(m,nrslab)=x2(2) 

xhc (m, nrslab - 1 )=gamma(m, nrs lab- l)*xj c(m, nrslab -1) 
xhs (m , nrslab - 1 ) -gamma (m , nr s lab  - 1 ) *xj s (m , nrs lab  - 1 ) 
do  n-nrslab-2 , 2 , -1 

num-s  sq  ( n+1 )  *  ( j  rho  1  ( m ,  n+1 )  -i-gamma  (m,n+l)*hrhol(m,n+l)) 
den“ssq(n)*(  jrho2(m,n)+gam]na(m,n)*hrho2(ni,n) ) 
x j  c (m , n ) -num*x j c(m,n+l)/den 
xhc  (in ,  n ) -ganuna  ( m ,  n )  *x j  c  ( m ,  n ) 
x  j  s  ( m ,  n )  -num*x  j  s  (in ,  n+1 )  /den 
xhs (m,n) -gamma (m,n)*xjs(m,n) 
end  do 

num-ssq(2)*(jrhol(m, 2)+gamma(m, 2)*hrhol(m, 2) ) 
den-ssq(l)*j rho2(m, 1) 
xjc(m, l)-num*xjc(m,2)/den 
xj  s (m , l)-num*xj  s (m, 2) /den 
end  do 
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go  to  2 
end  if 

if(nt  .eq.  2  .and.  nrslab  .gt.  3) then 
do  m-0 , mmax - 1 

a4(l,l)--ssq(l)*j  rho2 (m, 1) 
a4(l , 2)-ssq(2)*j rhol(m,2) 
a4(l , 3)-ssq(2)*hrhol(m, 2) 
a4(l,4)-0. 

a4(2 , l)--s(l)*jprho2 (m, 1) 
a4(2 , 2)-s (2)*jprhol(m, 2) 
a4(2 , 3)-s(2)*hprhol(m, 2) 
a4(2 , 4)-0 . 
a4 ( 3 , 1 ) ”0 . 

a4(3 , 2)-ssq(2)*j rho2(m, 2) 
a4(3 , 3)-ssq(2)*hrho2 (m, 2) 

a4 (3 ,4)-- ssq(3)*(gammal(m, 3)*jrhol(m, 3)+hrhol(m, 3) ) 
a4  ( 4 , 1 )  — 0  . 

a4 (4 , 2 ) -s ( 2 ) * j  prho2 (m , 2 ) 
a4(4 , 3)-s(2)*hprho2 (m,  2) 

a4(4,4)--s(3)*(gammal(m, 3)*jprhol(m, 3)+hprhol(m, 3) ) 

a4sav(l , l)-a4(l , 1) 

a4sav(l,2)-a4(l,2) 

a4sav(l, 3)-a4(l, 3) 

a4sav(l,4)-a4(l ,4) 

a4sav(2 , l)-a4(2 , 1) 

a4sav(2 , 2)-a/.  (2,2) 

a4sav(2 , 3)-a4(2 , 3) 

a4sav(2 ,4)-a4(2 ,4) 

a4sav(3 , l)-a4(3 ,1) 

a4sav(3 , 2)-a4(3 , 2) 

a4sav(3 , 3)-a4(3 , 3) 

a4sav(3 , 4)-a4(3 ,4) 

a4sav(4 , l)-a4(4, 1) 

a4sav(4 , 2)-a4(4 , 2) 

a4sav(4 , 3)-a4(4 , 3) 

a4sav(4 , 4)— a4(4 ,4) 

b4(l)~ssq(2)*ai"pc*hpr0(m)*j rhol (in,  2) 
b4(2)-s(2)*ampc*hpr0(m)*jprhol(m, 2) 
b4(3)-ssq(2)*ampc.*jpr0(m)*hrho2  Cm,  2) 
b4(4)-s(2)*ampc*jpr0(m)*hprho2 (m, 2) 
if(m  .eq.  0)then 
b4(l)-.5*b4(l) 
b4(2)-.5*b4(2) 
b4(3)-.5*b4(3) 
b4(4)-.5*b4(4) 
end  if 

call  clineq(a4,b4,x4,4,4,0,err) 
xjc(m, 1 ) — x4 ( 1 ) 
xjc(m, 2)-x4(2) 
xhc(m, 2)-x4(3) 
xhc(m, 3)-x4(4) 

b4(l)--ssq(2)*amps*m*gr0(m)*jrhol(m, 2) 
b4(2)--s(2)*amps*m*gr0(ni)*jprhol(m,  2) 
b4(3)--ssq(2)*amps*m*lr0(m)*hrho2 (m, 2) 
b4(4)--s(2)*amps*m*lr0(m)*hprho2(m, 2) 
call  clineq(a4sav,b4,x4,4,4,0,err) 
xjs(m,l)-x4(l) 
xjs(m,2)-x4(2) 
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xhs(m, 2)-x4(3) 
xhs(m, 3)-x4(4) 

xjc(m, 3)-gammal(m, 3)*xhc(m, 3) 
xjs(m, 3)-gammal(m, 3)*xhs(m, 3) 
do  n-4, nrslab-1 

num-ssq(n-l)*(gammal(m,n-l)*jrho2(m,n-i)+hrho2(m,n-l) ) 
den-ssq(n)*(gammal(m,n)*jrhol (m,n)+h- aol(m,n) ) 
xhc(m,n)-nuro*xhc(m,n-l) /den 
xjc(m,n)-gammal(m,n)*xhc(m,n) 
xhs(m,n)-num*xhs(m,n-l) /den 
xj s (m , n) -gammal (m , n) *xhs (m , n) 
end  do 

num-s  s  q(nrslab-l)*( gamma 1 ( m , nr  s 1 ab - 1 ) * j rho2 (m , nrslab- 1 ) + 

$  hrho2(m,nrsiab-l) ) 

den-ssq(nrslab)*hrhol(m, nrslab) 
xhc (m , nrslab) -num*xhc (m , nrslab - 1 ) /den 
xhs (m , nr slab)=num*xhs(m, nrslab-1) /den 
end  do 
go  to  2 
end  if 

if(nt  .eq.  nrslab-1  .and,  nrslab  ,gt.  3)then 
do  m-0 ,  nunax  - 1 

a4(l , l)-*-ssq(nrslab- 2)* (jrho2(m, nrslab- 2)+gamma(m, nrslab- 2)* 
$  hrho2(m,nrslab-2)) 

a4(l , 2 )“ssq( nrslab- l)*j rhol (m , nrslab-1) 
a4(1  ,  3~i - jsq '’nrslab -l)*hrhol (in .nrslab  - 1) 
a4(l,4)-0. 

a4(2,l)--s (nrslab- 2 )*( j prho2 (m , nrslab- 2)+gamma(m , nrslab- 2 ) * 

$  hprho2(ra,nrslab-2) ) 

a4 ( 2 , 2 ) -s ( nrslab - 1 ) * j  prhol (m , nrslab - 1 ) 
a4(2 , 3)-s(nrslab-l)*hprhol(m, nrslab-1) 
a4 ( 2 , 4 ) —0 . 
a4  ( 3 , 1 )  — 0 . 

a4(3 , 2)— ssq(nrslab- l)*j rho2(m,nrslab-l) 
a4(3 , 3 )-ssq (nrslab- l)*hrho2 (m, nrslab-1) 
a4(3 ,4)--ssq(nrslab)*hrhol(m,nrslab) 
a4 ( 4 , 1 ) -0 . 

a4 (4 , 2)-s (nrslab- l)*jprho2(m, nrslab-1) 

a4(4 , 3) -s (nr s lab -l)*hprho2(m, nrslab-1) 

a4 (4, 4) --s (nrslab) *hp rhol (m, nrslab) 

a4sav(l , l)-a4(l , 1) 

a4sav(l , 2)-a4(l , 2) 

a4sav(l , 3)-a4(l , 3) 

a4sav(l ,4)-a4(l ,4) 

a4sav(2 , l)-a4(2 , 1) 

a4sav(2 , 2)-a4(2 , 2) 

a4sav(2 , 3)-a4(2 , 3) 

a4sav(2 ,4)-a4(2 ,4) 

a4sav(3 , l)-a4(3 , 1) 

a4sav(3 , 2)-a4(3 , 2) 

a4sav(3 , 3)-a4(3 , 3) 

a4sav(3 ,4)-a4(3 ,4) 

a4sav(4 , l)-a4(4 , 1) 

a4sav(4 , 2)-a4(4 , 2 ) 

a4sav(4 , 3)*=a4(4 , 3) 

a4sav(4 ,4)-a4(4 ,4) 

b4(l)-ssq (nrslab - l)*ampc*hprO(m)*j rhol(m,nrslab- 1) 
b4(2)-s(nrslab-l)*ampc*hpr0(ra)*jprhol(m, nrslab- 1) 


b4(3)-ssq{nrslab- I)*ampc*jpr0(m)*hrho2(m,nrslab-1) 
b4(4)-s(nrslab-l)*ampc*jpr0(m)*hprho2(m,nrslab-l) 
if(m  .eq.  0)then 
b4(l)-.5*b4(l) 
b4(2)-.5*b4(2) 
b4(3)-.5*b4(3) 
b4(4)-.5*b4(4) 
end  if 

call  clineq(a4,b4,x4,4,410,err) 

xjc(m,nrslab-2)-x4(l) 

xj  c  (m ,  nr  s  lab  - 1 ) -x4 ( 2 ) 

xhc (m, nrslab -l)-x4(3) 

xhc(m,nrslab)-x4(4) 

b4(l)--ssq(nrslab-l)*amps*m*gr0(m)*j rhol(ra,nrslab-l) 

b4(2)--s(nrslab-l)*amps*m*gr0(m)*jprhol(m,nrslab-l) 

b4(3)--ssq(nrslab-l)*amps*m*lr0(m)*hrho2(m, nrslab -1) 

b4(4)--s (nrslab- I)*amps*m*lr0(m)*hprho2 (m, nr slab- 1) 

call  clineq(a4sav,b4,x4,4,4,0,err) 

xj  s (m , nr s lab  -  2 ) -x4 ( 1 ) 

xj  s (m , nrslab - 1 ) -x4 ( 2 ) 

xhs (m , nrslab - 1 ) -x4 ( 3 ) 

xhs (m,nrslab)-x4 (4) 

xhc(m,nrslab-2)-gamma(m,nrslab-2)*xjc(m,nrslab-2) 
xhs (m, nrslab -2)-gamma(m, nrslab- 2)*xj  s(ra, nrslab- 2) 
do  n— nrslab- 3 , 2 , -1 

num-ssq(n+l)*( jrhol(m,n+l)+gamma(m,n+l)*hrhol(m,n+l) ) 
den-ssq(n)*(  j  rho2(m,n)+gamma(m,n)*hrho2  (m,n) ) 
xj  c (m , n) -num*xj  c (m , n+1 ) /den 
xhc  (m ,  n)  -gamma (m ,  n)  *x  j  c  (m ,  n) 
xj s (m , n)-num*xj s (m , n+1 ) /den 
xhs (m , n) -gamma (m , n) *xj s (m , n) 
end  do 

num— ssq(2)*( j rhol(m,  2 ) -(gamma (m,  2)*hrhol (m ,  2) ) 
den-ssq(l)*j  rho2 (m, 1) 
xjc(m, l)-num*xjc(m, 2)/den 
xjs(m, l)-num*xjs(m, 2) /den 
end  do 
go  to  2 
end  if 


nt>4,nt  .ne.  1 , 2 .nrslab- 1 .nrslab 
do  m-0,mmax-l 

a4(l,l)--ssq(nt-l)*(j  rho2 (m, nt-l)+gamma(m, nt- l)*hrho2 (m.nt- 1) ) 
a4(i ,  2)-ssq(nt)*jrhol(iu,nt) 
a4(l , 3)-ssq(nt)*hrhol(m,nt) 
a4(l,4)-0. 

a4(2 , l)--s(nt-l)*(jprho2(m,nt-l)+gamina(m,nt-l)*hprho2(m,nt-l) ) 

a4(2 , 2)-s(nt)*jprhol(m,nt) 

a4(2 , 3)-s(nt)*hprhol(m,nt) 

a4 ( 2 , 4 ) -0 . 

a4(3 , l)-0 . 

a4(3 ,2)-ssq(nt)*jrho2(m,nt) 
a4(3 , 3)-ssq(nt)*hrho2(m,nt) 

a4(3 ,4)--ssq(nt+l)*(gammal(m,nt+l)*j rhol(m,nt+l)+hrhol(m,nt+l) ) 
a4(4,l}-0. 

a4(4, 2)-s(nt)^jprho2 (a.nt) 
a4(4, 3)-s(nt)*hprho2(m,nt) 

a4(4,4)--s(nt+l)*(gammal(m,nt+l)*jprhol(m,nt+l)+hprhol(m,nt+l) ) 
a4sav( 1 , 1 ) -a4 (1,1) 
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a4sav(l , 2)-a4(l , 2) 
a4sav(l , 3)-a4(l , 3) 
a4sav(l ,4)-a4(l ,4) 
a4sav(2 , l)-a4(2 , 1) 
a4sav(2 , 2)-a4(2 , 2) 
a4sav(2 , 3)-a4(2 , 3) 
a4sav(2 ,4)-a4(2 ,4) 
a4sav(3 , l)-a4(3 , 1) 
a4sav(3 , 2)-a4(3 , 2) 
a4sav(3 , 3)— a4(3 , 3) 
a4sav(3 ,4)— a4(3 ,4) 
a4sav(4 , l)-a4(4 , 1) 
a4sav(4 , 2)— a4(4 ,2) 
a4sav(4, 3)-a4(4, 3) 
a4sav(4,4)-a4(4,4) 

b4(l)-ssq(nt)*ampc*hprO(m)*jrhol(m,nt) 
b4(2)-s (nt)*ampc*hprO(m)*jprhol(m,nt) 
b4(3)-ssq(nt)*ampc*jpr0(m)*hrho2 (m ,nt) 
b4(4)-s(nt)*ampc*jpr0(m)*hprho2(m,nt) 
if(m  .eq.  0)then 
b4(l)-.5*b4(l) 
b4(2)-.5*b4(2) 
b4(3)“. 5*b4(3) 
b4(4)- . 5*b4(4) 
end  if 

call  clineq(a4,b41x4,4,4,0, err) 
xjc(m,nt-l) =x4 ( 1 ) 
xjc(m,nt)-x4(2) 
xhc(m,nt)-x4(3) 
xhc(m,nt+l)-x4(4) 

b4(l)--ssq(nt)*amps*m*grO(m)*jrhol(m,nt) 

b4(2)--s(nt)*amps*m*gr0(m)*jprhol(m,nt) 

b4(3)--ssq(nt)*amps*m*lr0(m)*hrho2(m,nt) 

b4(4)--s (nt)*amps*m*lr0(ra)*hprho2 (m,nt) 

call  clineq(a4sav,b4 , x4 ,4 ,4 ,0 , err) 

xjs(m.nt-l) -x4 ( 1 ) 

xjs(m,nt)-x4(2) 

xhs(m,nt)-x4(3) 

xhs ( m , nt+1 ) -x4 ( 4 ) 

xhc ( m , nt - 1 ) -gamma ( m , nt - 1 ) *x j  c ( m , nt - 1 ) 
xj c(m, nt+1 )-gammal(m, nt+1 )*xhc(m, nt+1) 
xhs (m , nt - 1 ) -gamma (m , nt - 1 ) *x j  s(m,nt-l) 
xj  s (m , nt+1 ) -garamal (m , nt+1 ) *xhs (m , nt+1 ) 
do  n-nt+2 ,nrslab-l 

num-ssq(n-l)*(gammal(m,n-l)*j rho2(m,n-l)+hrho2(m,n-l) ) 
den-ssq<n)*(gammal(m,n)*j rhol(m,n)+hrhol(m,n) ) 
xhc(m , n)-num*xhc (m , n- 1 ) /den 
xjc(m,n) -gammal (m,n)*xhc(m,n) 
xhs (m, n)-num*xhs(m, n- l)/den 
xj  s (m , n) -gammal (ra , n) *xhs (m , n) 
end  do 

num-ssq ( nr slab- 1)*( gammal (m,nrslab-l)*jrho2(m,nrs lab -1)+ 
hrho2(m,nrslab-l) ) 
den-  ssq(nrslab)*hrhol(ra,nrslab) 
xhc(m,nrslab) -num*xhc (m,nrslab-l)/den 
xhs(m,nrslab)-num*xhs(m,nrslab-l)/den 
do  n-r.t-2 ,2,-1 

n'UH-ssq(n+l)*(jrhol(m,n+l)+gamma(m,n+l)*hrhol(m. ”+1) ) 


m 


den-ssq(n)*(jrho2(m,n)+gamma(m,n)*hrho2(m,n)) 
xj  c (m , n)-num*xj  c (m , n+1 ) /den 
xhc (m , n) -gamma (m , n)*xj c (m , n) 
x j s ( m , n ) -num*x j s ( m , n+1 ) /den 
xhs (m , n) -gamma (m , n) *xj  s (m , n) 
end  do 

num-ssq(2)  +  ( j  rhol(m, 2)+gamma(m, 2)*hrhol(m, 2) ) 
den-ssq(l)*j rho2(m, 1) 
xj  c  (m ,  1 )  -num*j:j  c  (m ,  2 )  /den 
xjs(m,l)-num*xjs(m,2)/den 
end  do 
continue 
ntsav— nt 
return 
end 


subroutine  mlflds 
implicit  complex*16(a-h,o-z) 

complex*16  ks , im, jprhor , j rhor , j prhot , j rhot , jprO , j rO , IrO 
real*8  rmin, r , rsav, alpha, cosalp , s inalp , smin.xO ,xyint ,y0 , rO , rinc , 
$rmax , sine , smax , ctheta , rhorcv.pi , phi , s theta, rho , emphi , smphi , chi , 
$cchi , schi , ctmdl , ctmd2 , deltal , delta2 , stmdl , stmd2 , theta , scoord, 
$ezdb , ezang , hydb , hyang , hrdb , hrang , skip 
real*4  xplot ,yplotl ,yplot2 ,yplot3 , yplot4 ,yplot5 ,yplot6 
parameter (mmax-120) 
parameter (mxs lab-6 6) 
parameter (mxplo t-30 2) 

dimension  rho(mxslab) ,ks(mxslab) , s(mxslab) .ssq(mxslab) , 
$hprhor(0:mmax) ,hrhor (0 :mmax) , jprhor (0:mmax) , jrhor(0 :mmax) , 
$xjs(0:mmax,mxslab) ,xjc(0:mmax,mxslab) , xhs (0 :mmax,mxslab) , 
$xhc(0:mmax,mxslab) , j prhot (0 :mmax) , jrhot(0:mmax) ,hprhot(0 :mmax) , 
$hrhot(0 :mmax) , jr0(0:mmax) , jpr0(0:mmax) ,hr0(0:mmax) ,hpr0(0 :mmax) , 
$bjr0(0:mmax) ,bjpr0(0:mmax) ,bhr0(0 :mmax) ,bhpr0(0 :mmax) , 
$bgr0(0:mmax) ,blr0(0 :mmax) , 

$gr0(0:mmax) ,lr0(0:mmax) .xplot(mxplot) .yplotl(mxplot) , 
$yplot2(mxplot) ,yplot3(mxplot) ,yplot4(mxplot) ,yplot5(mxplot) , 
$yplot6 (mxplot) 

common/one/ampc , amps , ks .nrslab , rO , rho , s , ssq 

common/ three/alpha , amppre , ampprs , cl , c2 , cosalp , deltal , delta2 , 

$ if lag, rinc , rmax, rmin, s Inalp , sine , smax, smin.xyint , theta.xO ,y0 , skip 
common/ f our/xhc , xhs , xj  c , xj  s 
common/f ive/grO ,hpr0 ,hr0 , jprO , jrO , IrO 
data  pi/3.1415926536/ 
data  im/(0. , 1 . )/ 
r-rmin 

if (r  .eq.  0.)r-l.d-4 
rsav-0 . 
rtsav-nt 
scoord-smin 
s theta-sin (theta) 
c theta-cos ( theta) 
if(iflag  .eq.  0)then 

numpts-(rmax-rmin) /rinc+1 
write (40)  x0,y0,numpts 
else 

numpts-(smax-smin) /sinc+1 
write(40)  r.numpts 
end  if 

do  num— l.numpts 

rhorev-sqr t ( r*r+  r0*r0 -  2 . *r*r0*c  the  ta) 
if(rhorcv  .It.  l.)rhorcv-l. 
if (rO*rO+rhorcv*rhorcv-r*r  .eq.  0  )then 
phi-pi/2 . 

if(stheta  .It.  0.)phi--phi 
else 

phi-atan2(2 . *r*r0*s theta , rO*rO+rhorcv*rhorcv-r*r) 
end  if 
nr-nrslab 
do  n-l,nrslab-l 

if(rhorcv  . le .  rho (n)) then 
nr-n 
go  to  1 


end  if 
end  do 

1  continue 
nt-nrslab 

do  n-l.nrslab-l 

if(rO  .le.  rho(n))then 
nt-n 
go  to  2 
end  if 
end  do 

2  continue 

if(abs(r-rO)  .It.  150. )go  to  3 
if (abs(rhorcv-rO)  .It.  skip)go  to  3 
argr0-ks(nr)*r0 

call  hj  func ( argrO , b j  rO , b j  prO , bhrO , bhprO ) 
do  m— 0 ,  mraax  - 1 

blr0(m)-bjr0(m) /argrO 
bgr0(m)-bhr0(m) /argrO 
end  do 

if(r  .ne.  rsav)then 
arg-ks(nrslab)*r 
call  cbesjy (arg , 0 ,bj ,by , 0 , 0) 
hO-b  j - im*by 

call  cbesjy (arg, 1 ,bj ,by , 0 , 0) 
hl-bj - im*by 

call  cbesjy (arg, 2 ,bj ,by , 0 , 0) 
h2— b  j - im*by 
hpl— 0 . 5* (hO -h2 ) 
ezamb- ssq (nrs lab ) *ampprc*hl 
hy amb- - im* s ( nr s 1 ab ) *ampp r c*hp 1 
hramb— im*s(nrslab)*ampprs*hl/arg 
end  if 

if(nt  .ne.  ntsav  .or.  r  .ne.  rsav)then 
arg— ks(nt)*r 

call  cbesjy(arg,0,bj , by, 0,0) 
hO-b j - im*by 

call  cbesjy (arg, 1 ,bj ,by , 0 , 0) 
hl-bj - im*by 

call  cbesjy (arg, 2 ,bj , by , 0 , 0) 
h2-bj - im*by 
hpl-0 . 5*(h0-h2) 
eztst-ssq(nt)*ampprc*hl 
hy ts  t- - im*s ( nt ) *ampprc*hp 1 
hrtst-im*s (nt)*ampprs*hl/arg 
end  if 
ntsav-nt 

c  if (nr  .eq.  nt)then 

arg-ks(nr)*r 

call  cbesjy (arg, 0,bj , by, 0,0) 
hO-b j - im*by 

call  cbesjy(arg, 1 ,bj , by, 0,0) 
hl-bj - im*by 

call  cbesjy(arg, 2 ,bj , by, 0,0) 
h2-bj - im*by 
hpl-0. 5*(h0-h2) 
ezprmy-ssq(nr)*ampprc*hl 
hyprmy-- im*s (nr ) *ampprc*hpl 
hrpnDy-im*s(nr)*ampprs*hl/arg 
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end  if 

ar g-ks (nt ) *rhorcv 
ez sum-0 . 
hosum-0 . 
hrsum-0 . 

call  hj  func (arg , j  rhot , j  prhot ,hrhot , hprhot) 
do  m-mmax -1,0, - 1 
cmphi-cos (m*phi) 
smphi-sin(m*phi) 
if(rhorcv  .It,  r0)then 
if(m  .eq.  0)then 

ezsum-ezsum-ampc* .  5*ssq(iit)*hpr0(m)J'j  rhot (m)*cmphi 
hosum-hosum+ampc* . 5*im*s (nt ) *hpr0 (m) *j prhot (m) *cmphi 
else 

ezsum-ezsum-ampc*s£q(nt)*hprO(m)*j rhot(m)*cmphi 
$  +amps*ssq(nt)*m*grO(m)*j  rhot(m)*smphi 

hosum-hosum+ampc*im*s (nt ) *hpr0 (m) * j  prhot (m) *cmphi 
$  - amps*im*s(nt)*m*grO(m)*j prhot (m)*smphi 

hrsum=hrsum+ampc*im*s(nt)*m*hprO(m)*j rhot (m)*smphi/arg 
$  +amps*im*s (nt)*m*m*grO(m)*j  rhot(m)*cmphi/arg 

end  if 
else 

if(m  .eq.  0)  then 

ezsum-ezsum-ampc* . 5*ssq(nt)*jpr0(m)*hrhot(m)*cmphi 
hosum-hosum+ampc* . 5*im*s (nt ) * j  prO (m) *hprho t (m) *cmphi 
else 

ezsum-ezsum-ampc*ssq(nt)*jprO(m)*hrhot (m)*cmphi 
$  +amps*ssq(nt)*m*lrO(m)*hrhot(m)*smphi 

hosum-hosum+ampc*im*s (nt ) *j  prO (m) *hprhot (m) *cmphi 
$  -amps*im*s (nt)*m*lrO(m)*hprhot(ra)*smphi 

hrsum-hrsum+ampc*im*s (nt ) *m*j  prO (m) *hrhot (m) *smphi/arg 
$  *-amps*im*s  (nt )  *m*m*lr0  (m)  *hrhot  (m)  *cmphi/arg 

end  if 
end  if 
end  do 

chi-pi- theta-phi 

cchi-cos(chi) 

schi-sin(chi) 

hysum-hosum*cchi  -hrsum*sc.hi 
hrsum-hosum*schi+hrsum*cchi 
write (44,*) 'nt-' ,nt, 'nr-' ,nr 

write (44 ,*)' chi-' , chi , ' theta-' , theta , 'phi- ' , phi 
write (44,*) 'r-' ,r, 'r0=' ,r0, 'rhorcv-' .rhorcv 
write (44,*) 'x0-' ,x0, 'y0-' ,y0 

if (zabs (eztst-ezsum)  .gt.  .005  .or.  zabs(hytst-hysum)  -gt.  .005 
$  .or.  zabs (hrtst-hrsum)  .gt.  .005)then 
write(44 ,*)' convergence  warning' 
write(44,*)'eztst=' ,eztst, 'ezps-' ,ezsum 
write(44,*) 'hytst-' ,hytst, 'hyps-' .hysurn 
write (44 , *) 'hrtst-' , hrtst , ' hrps- ' , hr sum 
end  if 

arg-ks (nr)*rhorcv 

call  hj  func (arg , j  rhor , jprhor , hrhor ,hprhor) 
ezs-0 . 


write (44 ,*) 'm-' ,m, 'nr-' ,nr, 'xhc-' ,xhc(m,nr) , 'xhs-' ,xhs(m,nr) 
cmphi-cos (m*phi ) 
smphi-s in (m*ph 1 ) 
if(nt  .eq.  nr) then 
if (nr  .eq.  l)then 

temp j -(xj  c (m , nr ) *cmphi+xj  s (m , nr ) *smphi ) 
tempi j -m* ( - xj  c (m , nr ) *smphi+xj  s (m , nr ) *cmphi ) 
ezs-ezs+ssq(nr)*tempj*jrhor(m) 
hphis-hphis - im*s (nr ) *temp j  *j  prhor (m) 
hrhos-hrhos+im*s(nr)*templj*j rhor(m) /arg 

else 

if(nr  .eq.  nrslab)then 

temph-(xhc (m , nr) *cmphi+xhs (m , nr ) *smphi ) 
templh-m* ( -xhc (m , nr ) *smphi+xhs (m , nr ) *cmphi ) 
ezs-ezs+ssq(nr)*temph*hrhor (m) 
hphis-hphis- im*s (nr) *temph*hp rhor (m) 
hrhos-hrhos+im*s(nr)*templh*hrhor(m) /arg 

else 

temp j - (xj  c (m , nr ) *cmphi+xj  s ( m , nr ) *smphi ) 
temph- (xhc (m , nr ) *cmphi+xhs (m , nr ) *smphi ) 
tempi j -m* ( - x j  c ( m , nr ) *  smph i+x j  s ( m , nr ) *cmph i ) 
templh-m* ( -xhc(m,nr)*smphi+xhs(m,nr)*cmphi) 
ezs— ezs+ssq (nr ) * ( temp j  *j  rhor (m) +temph*hrhor (m) ) 
hphis-hphis - im*s (nr)* 

(tempj *j prhor (m)+temph*hprhor(m) ) 
hrhos-hrhos+im*s(nr)* 

( tempi j*j rhor(m)+templh*hrhor(m) ) /arg 

end  if 
end  if 

go  to  10 

end  if 

if(nr  .eq.  l)then 

if (rhorcv  .It.  r0)then 
if(m  .eq.  0)then 

tempj - ( x j  c ( m , nr ) *cmph i+x j  s ( m , nr ) *smph i ) 
temp I j -m* ( - x j  c ( m , nr ) *smph i+x j  s ( m , nr ) *cmph i ) 
ezs-ezs+ssq(nr)*tempj*jrhor(m) 

+ampc*. 5*ssq(nr)*bhpr0(m)*j rhor(m)*cmphi 
hphis-hphis-im*s(nr)*tempj*jprhor(m) 

- ampc* . 5* im*s ( nr ) *bhpr 0 (m) * j  prhor (m) *cmphi 
hrhos-hrhos+im*s ( nr ) *templ j  *j rhor (ra) /arg 
else 

tempj -(xjc(m,nr) *cmphi+x j  s ( m , nr ) *  smph i ) 
tempi j -m* (-xjc(m,nr) *  smph i+x j  s ( m , nr ) *cmph i ) 
ezs-ezs+ssq(nr)*tempj *j rhor(m) 

+ampc*ssq(nr)*bhprO(m) 'jrhor(m)*cmphi 
-amps*ssq(nr)*m*bgrO(m)*jrhor (m)*smphi 

hphis-hphis - im*s (nr ) * temp j  * j  prhor (m) 

- ampc*im*s (nr ) *bhpr0 (m) * j  prhor (m) *cmphi 
+amps*im*s (nr ) *m*bgrO (m) *j  prhor (m) *smphi 
hrhos-hrhos+im*s(nr)*templj *jrhor(m)/arg 

- ampc*im*s (nr ) *m*bhpr0 (m) *j  rhor (m) *smphi/arg 
-amps*im*s(nr)*m*m*bgrO(m)*jrhor(m)*cmphi/arg 

end  if 
else 

if(m  ,eq.  0)then 

tempj-(xjc(m,nr) *cmph i+xjs(m,nr)*s  mph i ) 
tempi j -m* ( -xj  c (m , nr ) *smphi+xj  s (m , nr ) *cmphi ) 


ezs-ezs+ssq(nr)*tempj*j  rhor(m) 

+ampc*. 5*ssq(nr)*bjpr0(m)*hrhor (m)*cmphi 
hphis-hphis - im*s (nr ) *temp j  *j prhor (m) 

-ampc* . 5*ira*s(nr)*bjpr0(m)*hprhor (m)*cmphi 
hrhos-hrhos+im*s (nr)* tempi j*j rhor (m) /arg 
else 

tempj-(xj  c(m,nr)*cmphi+xj  s (ra, nr)*smphi) 
templj-m*( -xjc(m,nr)*smphi+xj  s(m,nr)*cmphi) 
ezs-ezs+ssq(nr)*tempj*jrhor (m) 

+ampc*ssq(nr)*bjprO(m)*hrhor  (m)*cmp'ni 
-amps*ssq(nr)*m*blrO(m)*hrhor  (m)*smphi 
hphis-hphis - im*s (nr ) *temp j  *j prhor (m) 

- ampc*im*s (nr ) *b j  prO (m) *hprhor (m) *cmphi 
+amp  s* im*  s ( nr ) *m*b 1 r 0 ( m ) *hp  rhor ( m ) *smph i 
hrhos=hrhos+im*s (nr)*templj*j rhor (m) /arg 

-ampc*im*s(nr)*m*bjprO(m)*hrhor (m)*smphi/arg 
-amps*ira*s (nr)*m*m*blrO(m)*hrhor (m)*cmphi/arg 

end  if 
end  if 
se 

if (nr  .eq.  nrslab)then 
if(rhorcv  .It.  rO)then 
if(m  .eq.  0)then 

teraph- ( xhc ( m , nr ) *cmphi+xhs (m , nr ) *smph i ) 
templh-m*(  -xhc  (m,  nr)*smphi+xhs(ni,nr)*cmphi) 
ezs-ezs+ssq(nr)*temph*hrhor (m) 

+ampc*. 5*ssq(nr)*bhpr0(m)*j rhor (ra)*cmphi 
hphis-hphis - im^  j (nr)*temph*hprhor (m) 

-ampc*. 5*im*s(nr)*bhpr0(m)*jprhor(m)*cmphi 
hrhos-hrhos+im*s(nr)*templh*hrhor(m)/arg 
else 

temph-(xhc(m,nr)*cmphi+xhs(m,nr)*smphi) 
templh-m*( -xhc (m,nr)*smphi+xhs (m,nr)*cmphi) 
ezs-ezs+ssq(nr)*temph*hrhor (m) 

+ampc*ssq (nr ) *bhprO (m) *j rhor (m) *cmphi 
-amps*ssq(nr)*m*bgrO(m)*j rhor (m)*smphi 
hphis-hphis - im*s (nr )*temph*hprhor (m) 

-ampc*im*s ( nr )*bhprO(m)*j prhor (m)*cmphi 
+amps*im*s (nr )*m*bgrO(m)*j prhor (m)*smphi 
hrhos-hrhos+im*s(nr)*templh*hrhor (m)/arg 

- ampc* im*s (nr ) *m*bhprO (m) * j  rhor (m) *smphi /arg 
- amps*im*s ( nr ) *m*m*bgrO (m) * j rhor (m) *cmphi/arg 

end  if 
else 

if(m  .eq.  0)then 

temph-(xhc(m, nr)*cmphi+xhs (m,nr)*smphi) 
templh-ra* ( -xhc (m , nr ) *smphi+xhs (m , nr ) *cmphi ) 
ezs-ezs+ssq(nr)*femph*hrhor (m) 

+ampc* . 5*ssq ( nr ) *bj  prO (m) *hrhor (m) *cmphi 
hphis-hphis - im*s (nr)*temph*hprhor (m) 

- ampc* . 5*im*s (nr ) *b j  prO (m) *hprhor (m) *cmphi 
hrhos-hrhos+im*s (nr )*templh*hrhor (m) /arg 
else 

t emph- ( xhc ( m , nr ) +cmph i+xhs (m,nr) *smph i ) 
templh-m* ( - xhc (m , nr ) *smph i+xhs (m , nr )*cmphi ) 
ezs-ezs+ssq(nr)*temph*hrhor (m) 

+ampc*ssq(nr)*bjprO(m)*hrhor(m)*cmphi 
-amps*ssq(nr )*m*blrO(m)*hrhor (m)*smphi 
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hph i s-hph is-im*s(nr) *  t emph*hp rho r ( m ) 

-ampc*im*s(nr)*bj  prO(m)*hprhor (m)*cmphi 
+amps*im*s(nr)*m*blrO(m)*hprhor(m)*smphi 
hrhos-hrhos+im*s(nr)*templh*hrhor (m)/arg 

- ampc*im*s (nr ) *m*bj prO (m) *hrhor (m) *smphi/arg 
-amps*im*s  (nr)*m*m*blrO(m)*hrhor (m)*cmphi/arg 

end  if 
end  if 
else 

if(rhorcv  .It.  rO)then 

if(m  .eq.  0)then 

tempj-(xjc(m,nr)*cmphi+xj  s(m,nr)*smphi) 
temph-(xhc (m , nr ) *cmph i+xhs (m , nr )*smphi ) 
templj-m*(-xjc(m,nr)*smphi+xjs(m,nr)*cmphi) 
temp lh-m* ( - xhc ( m , nr ) *  smph i+xhs (m , nr ) *cmph i ) 
ezs-ezs+ssq(nr)*(tempj *j  rhor (m)+temph*hrhor (m) ) 
+ampc* . 5*ssq(nr)*bhpr0(m)*j  rhor(m)*cmphi 
hphis-hphis - im*s (nr)* 

( temp j  * j  prhor (m) +temph*hprhor (m) ) 

-ampc* . 5*im*s (nr )*bhprO(m)*j prhor (m)*cmphi 
hrhos=hrhos+im*s (nr)* 

( templj*j  rhor (m)+ temp lh*hrhor(m) ) /arg 

else 

t emp  j  - ( x j  c ( m , nr )  *cmph  i+x  j  s  ( m ,  nr )  *smph  i ) 
temph- (xhc(m,nr) *cmphi+xhs ( m , nr ) *smph i ) 
temp  1  j  -m*  ( -  x  j  c  ( m ,  nr )  *  smph  i+x  j  s  (m ,  nr )  *cmph  i ) 
temp lh=m* ( - xhc (m , nr ) *smphi+xhs (m , nr ) *cmphi ) 
ezs-ezs+ssq(nr)*(tempj *j  rhor (m)+temph*hrhor (m) ) 
+ampc*ssq (nr ) *bhprO (m) *j rhor (m)*cmphi 
- amps*ssq (nr ) *m*bgrO (m) *j  rhor (m) *smphi 
hphis-hphis - im*s (nr)* 

( temp j  *j prhor (m)+temph*hp rhor (m) ) 

-ampc*im*s (nr) *bhprO(m)*j prhor (m)*cmphi 
+amps*im*s (nr ) *m*bgrO (m) *j  prhor (m) *smphi 
hrhos-hrhos+im*s(nr)* 

(templj  *j  rhor (m)+templh*hrhor(m) )/arg 
- ampc*im*s (nr ) *m*bhprO (m) *j  rhor (m) *smphi/arg 
-amps*im*s (nr)*m*m*bgrO(m)*j  rhor (m)*cmphi/arg 

end  if 
else 

if(m  .eq.  0)then 

tempj-(xjc(m,nr)*cmphi+xj  s(m,nr)*smphi) 
temph- ( xhc ( m , nr ) *cmphi+xh  s ( m , nr ) *smph i ) 
temp 1 j -m* ( - x j  c ( m , nr ) *smph i+x j  s ( m , nr ) *cmph i ) 
temp lh-m* ( -xhc (m , nr )*smph i+xhs (m , nr ) *cmphi) 
ezs-ezs+ssq(nr)*(tempj*j rhor (m)+temph*hrhor(m) ) 
+ampc* . 5*ssq(nr)*bjpr0(m)*hrhor (m)*cmphi 
hphis-hphis - im*s(nr)* 

( temp j  * j  prhor (m) +temph*hprhor (m) ) 

-ampc*. 5*im*s (nr)*bjprO(m)*hprhor(m)*cmphi 
hrhos-hrhos+im*s (nr)* 

(templj  *j rhor (m)+templh*hrhor (m) ) /arg 

else 

t emp j -(xjc(m,nr) *cmph i+x j  s ( m , nr ) *smph i ) 
temph- (xhc (m , nr ) *cmphi+xhs (m , nr )*smphi ) 
templj -m* ( - x j  c ( m , nr ) *smph i +x j  s ( m , nr ) *cmph i ) 
temp lh-m* ( -xhc (m , nr ) *smphi+xhs (m , nr )*cmphi ) 
ezs-ezs+ssq(nr)*( tempj *j rhor (m)+temph*hrhor (m) ) 
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+ampc*ssq(nr)*bjprO(m)*hrhor (m)*cmphi 
-amps*ssq(nr )*m*blrO(m)*hrhor(m)*smphi 
hphis-hphis - im*s (nr )* 

(tempj*jprhor (m)+temph*hprhor(m) ) 

- ampc*im*s ( nr ) *b j prO (m) *hprhor (m) *cmphi 
+amps*im*s (nr )*m*blrO(m)*hprhor (m)*smphi 
hrhos-hrhos+im*s (nr ) * 

( tempi j *j rhor (m) +templh*hrhor (m) ) /arg 
'  -ampc*im*s (nr)*m*bjprO(m)*hrhor (m)*smphi/arg 

i  -amps*im*s(nr)*m*ni*blrO(in)*hrhor  (m)*cmphi/arg 

end  if 
end  if 
end  if 
end  if 
continue 
end  do 

chi-pi- theta -phi 
cchi=cos (chi) 
schi-sin(chi) 
ez-ezprmy+ezs 

hy=hyp  rmy+hph i s  *c  ch i - hrho  s  *s  ch i 
hr=hrprmy+hphis*schi+hrhos*cchi 
write(44 ,*) ' scoord=' .scoord 
write (44 , *) ' ezamb- ' , ezamb 
if (zabs (ampprc)  .ne.  0.)then 
qez-zlog(ez/ezamb) 
qhy-z 1 og ( hy /hy amb ) 
end  if 

if (zabs(ampprs)  .ne.  0 . )qhr-zlog(hr/hramb) 

ezdb-8 . 6858896*qez 

hydb-8 . 6858896*qhy 

hrdb-8 . 6858896*qhr 

ezang—- im*qez*57 . 295779 

hyang-- ira*qhy*57 . 295779 

hrang=- im*qhr*57 .295779 

if(ezang  .gt.  180 . )ezang-ezang- 360 . 

if (ezang  .It.  - 180 . )ezang=ezang+360 . 

ezang— -ezang 

if(hyang  .gt.  180. )hyang»hyang-360. 
if(hyang  .It.  -180 . )hyang=hyang+360 . 
hyang--hyang 

if(hrang  .gt.  180 . )hrang-hrang- 360 . 
if(hrang  .It.  -180 . )hrang-hrang+360. 
hrang--hrang 
yplotl (num)-sngl (ezdb) 
yplot2 (nura) -sngl (hydb) 
yplot3 (num)-sngl (hrdb) 
yplot4(num) -sngl (ezang) 
yplot5 (num)-sngl (hyang) 
yplot6(num)=sngl(hrang) 
if(iflag  ,eq.  0)then 
xplot(num)-sngl (r) 
else 

xplot(num) -sngl (scoord) 
end  if 

write(44,*) ' ez/ezu(db)-’ ,ezdb 
write  (44  ,*) '  ez/''zu(deg)-'  ,  ezang 
write (44,*) ' hy/hyu ( db ) - ' ,hydb 
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write(44,*) 'hy/hyu(deg)-' ,hyang 
write(44,*) *hr/hru(db)-' ,hrdb 
write(44,*) 'hr/hru(deg)-' .hrang 

write(40)  xplot(num) ,yplotl(num) ,yplot2(num) ,yplot3(num) , 
$  yplot4(num) ,yplot5(num) ,yplot6 (num) 

continue 

if (if lag  .eq.  0)then 
r-r+rinc 

if(r  .eq.  0.)r-l.d-4 
if(r  .gt.  rmax+.01)go  to  5 
else 

if (alpha  .eq.  0.)then 
xO-xO+sinc 
yO-xyint 

scoord— scoord+s inc 
else 

xO-xO+sinc*cosalp 
yO-yO+sinc*sinalp 
scoord-scoor  d+s inc 
end  if 

if(scoord  .gt.  smax+.01)go  to  5 
rO-sqrt (xO**2+yO**2 ) 
if (rO  .It.  1 . ) then 
if(xO  .It.  0.)then 
xO--l . 

rO-sqrt (xO**2+yO**2) 
else 
xO-1 . 

rO— sqrt(xO**2+yO**2) 
end  if 

write(44,*) 'r0“' ,r0 
go  to  4 
end  if 

if(xO  .eq.  0.)then 
theta-pi/2 . 

if(yO  .It.  0. ) theta-- theta 
else 

theta-atan2 (yO , xO) 
end  if 

ctmdl-cos(theta-deltal) 
ctmd2-cos(theta-delta2) 
stmdl-sin( theta- deltal) 
s  tmd2-s in ( the  ta - de 1 ta2 ) 
s  the  ta-s in ( theta) 
ctheta-cos( theta) 
ampc-2 .*ctmdl*cl+2 .*ctmd2*c2 
amps-2 . *stmdl*cl+2 . *stmd2*c2 
call  mlcoef 
end  if 
end  do 
return 
end 
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SUBROUTINE  CBESJY(Z , K , BJ , BY , KIND , NPRINT) 

IMPLICIT  C0MPLEX*16  (A-H.O-Z) 

REAL* 8  DGMSUM , DG1 , DG2 , PI , RRT , EULER , KFAC , KM1FAC 

data  PI/3. 14159265358979D0/, EULER/0. 577215664901533D0/ 

C 

IF( ZABS (Z)  .NE.  0 . DO)  GO  TO  7 
BJ-(0. 000,0. 0D0) 

IF(K  .EQ.  0)  BJ-(1 . 0D0 , 0 . 0D0) 

IF(KIND  .NE.  I  .AND.  K  . EQ .  0)  PRINT  400 
400  FORMAT ( 1H0, '***  Y  NOT  CALCULATED  FOR  ARGUMENT  OF  MAGNITUDE  0'//) 
RETURN 

7  IF(ZABS(Z)  .LT.  13. DO)  GO  TO  10 
C 

ASSYMPTOTIC  EXPANSION 

RHO-8 . *Z 
MU~4*K**2 

RT-ZSQRT(2./(PI*Z)) 

RRT-RT 

IF(RRT  .LT.  0.0D0)  RT- - RT 
P-0. 

C  DO  LOOP  FOR  CALCULATING  P 
DO  1  N-1,30 
M-N-l 
MM-2*M 

IF(N  .EQ.  1)  GO  TO  2 

TERM— ( -l)*(MU-( 2*MM- 3 )**2 )*(MU- (2*MM-l)**2)/(MM*(MM-l)*RHO**2)* 

$  TEPM 

IF (NPRINT  .EQ.  1)  PRINT  100 .TERM 
100  FORMAT ( '  TERM-' ,2030.15) 

P-P+TERM 

IF(ZABS(TERM)  .GE.  ZABS (TERMS)  .OR.  ZABS (TERM)  .LE.  l.D-17) 

$  GO  TO  3 
TERMS-TERM 
GO  TO  1 

2  TERM-1 . DO 
P-TERM 

IF(NPRINT  .EQ.  1)  PRINT  100, TERM 
TERMS-TERM 
1  CONTINUE 

3  CONTINUE 

IF(NPRINT  .EQ.  1)  PRINT  200, P 
200  FORMAT (1H0 , ' P-' , 2D30 .15) 

Q-0. 

C 

C  DO  LOOP  FOR  CALCULATING  Q 
DO  4  N-1,30 
M-N-l 
MM-2*M 
MMM-2*M+1 

IF(N  .EQ.  1)  GO  TO  5 

TERM-( -1)*(MU- (2*MM-1)**2)*(MU- (2*MM+1)**2)/(MMM*(MMM- l)*RMO**2)* 
$  TERM 

IF(NPRINT  .EQ.  1)  PRINT  100, TERM 
Q-Q+TERM 

IF(ZABS(TERM)  .GE.  ZABS (TERMS)  .OR.  ZABS (TERM)  .LE.  l.D-17) 

$  GO  TO  6 


TERMS-TERM 
GO  TO  4 

5  TERM-(MU-l)/RHO 
Q-TERM 

IF(NPRINT  .EQ.  1)  PRINT  100, TERM 
TERMS-TERM 
4  CONTINUE 

6  CONTINUE 

IF(NPRINT  .EQ.  1)  PRINT  300, Q 
300  FORMAT ( 1H0, 'Q-' ,2D30. 15) 

BJ— RT*ZCOS (Z-K*PI/2 . -PI/4 .) *P-RT*ZSIN(Z-K*PI/2 . -PI/4 .) *Q 
IF(KIND  .EQ.  1)  GO  TO  8 

BY-RT*ZSIN(Z-K*PI/2 . -PI/4. )*P+RT*ZCOS(Z-K*PI/2 . -PI/4. )*Q 
8  RETURN 

C 
C 

C  POWER  SERIES  EXPANSION 
C 

10  NTERMS-35 

KFAC-1 

IF(K  . LE.  1)  GO  TO  30 
DO  20  N-2 ,K 
20  KFAC-KFAC*N 

30  TERM— (Z/2 . )**K/KFAC 

BJ-TERM 

IF(KIND  .EQ.  1)  GO  TO  91 
DG1-0 . 

DG2-0 . 

IF(K  .EQ.  0)  GO  TO  80 
DO  60  N-l.K 
60  DG2-DG2+1./N 

80  TSUM3— TERM*DG2 

SUMT3— TSUM3 

91  DO  40  M-l , NTERMS 

TERM-TERM* ( Z/2 . ) **2* ( - 1 ) / ( (K+M) *M) 

BJ-BJ+TERM 

IF(KIND  .EQ.  1)  GO  TO  92 
DG1-DG1+1 . DO/M 
DG2— DG2+1 . D0/(M+K) 

DGMSUM-DG 1+DG2 
TSUM3— TERM*DGMSUM 
SUMT3-SUMT3+TSUM3 

92  IF(ZABS(TERM)  .LE.  l.D-17)  GO  TO  50 

40  CONTINUE 

50  IF(KIND  .EQ.  1)  GO  TO  93 

TERM3-SUMT3/PI 

TERMl-(2 ,/PI)*BJ*(EULER+ZLOG(Z/2 . ) ) 

SUMT2— (0 . ,0.) 

IF(K  .EQ.  0)  GO  TO  120 
KM1FAC-KFAC/K 
TSUM2-KM1FAC/ ( (Z/2 . )**K) 

SUMT2-TSUM2 

IF(K  .EQ.  1)  GO  TO  120 
KM1-K-1 
DO  130  M-l.KMl 
KMM-K-M 

TSUM2-TSUM2/(KMM*M)*(Z/2 . )**2 
130  SUMT2-SUMT2+TSUM2 


TERM2--SUMT2/PI 

BY-TERM1+TERM2+TERM3 

RETURN 

END 


SUBROUTINE  HJFUNC(ARG , BJ , BJP ,H2 ,H2P) 


C 

C  SPECIAL  PURPOSE  ROUTINE  FOR  CALCULATING  BESSEL  FUNCTIONS ( JM, VM) , 

C  HANKEL  FUNCTIONS (H2)  AND  THEIR  DERIVATIVES  FOR  COMPLEX  ARGUMENT. 

C  REQUIRES  FOR  STARTING  CONDITIONS  AN  INDEPENDENT  CALCULATION  OF 

C  JO, YO, J1  AND  Yl.  THE  LATTER  ARE  OBTAINED  FROM  SUBROUTINE  CBESJY. 

C 

IMPLICIT  COMPLEX*16 (A-H , 0 - Z) 

COMPLEX*16  IM 
REAL* 8  PI 

PARAMETER  (MMAX-51) 

PARAMETER  (MPLUS-MMAX+30) 

DIMENSION  BJ (0 :MMAX) , BJP(0 :MMAX) ,H2(0:MMAX) ,H2P(0:MMAX) ,BY(0:MMAX) 
DATA  IM/(0. 0,1.0)/ 

DATA  PI/3 . 14159265358979D0/ 

C 

C  CALCULATE  BY 

M-0 

CALL  CBESJY(ARG ,M, BJEXAC , BY(M) ,0,0) 

M-l 

CALL  CBESJY(ARG,M,BJTEMP,BY(M) ,0,0) 

DO  100  M-2.MMAX 

BY(M)— 2 . 0*(M- 1) /ARG*BY(M- 1) -BY(M-2) 

100  CONTINUE 
C 

C  CALCULATE  BJ 

BJP2-0.0 
BJP1-1.0 

DO  150  M-MPLUS-2 .MMAX-l , -1 
BJTEMP-2 . 0* (M+l ) /ARG*BJ  PI - B J  P2 
BJP2-BJP1 
BJP1-BJTEMP 

IF(ZABS (BJPl)  .GT.  1.D20)  THEN 
BJP2— BJP2/BJP1 
BJPl— 1 . 0 
END  IF 

150  CONTINUE 
C 

BJ (MMAX)— BJP2 
BJ (MMAX- 1)— BJPl 
C 

DO  200  M-MMAX-2 ,0,-1 

BJ (M) -2 . 0* (M+l ) /ARG*BJ (M+l ) - BJ (M+2 ) 

IF(ZABS(BJ (M) )  .GT.  1.0D20)  THEN 
BJSAVE-BJ (M) 

DO  199  MM-M.MMAX 

199  BJ (MM)-BJ (MM) /BJSAVE 
END  IF 

200  CONTINUE 

TERM— BJEXAC/BJ (0) 

DO  300  M-0, MMAX 
BJ (M)-TERM*BJ (M) 

300  CONTINUE 
C 

C  CALCULATE  H2 

DO  350  M-0, MMAX 
H2 (M)-BJ (M) - IM*BY(M) 

350  CONTINUE 
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CALCULATE  BJP  AMD  H2P 
BJP(O)- •  B  J  ( 1 ) 

H2P(0)  — H2(l) 

DO  400  M-l.MMAX-l 

BJP(M)— 0 . 5*(BJ(M- 1) -BJ (M+l) ) 

H2P(M)-0 . 5*(H2(M-1) -H2 (M+l) ) 

CONTINUE 


WEXACT-2 . 0/  (  PI*ARG  ) 

WMEQ0=B J ( 1 ) *BY ( 0 ) - BJ ( 0 ) *BY ( 1 ) 

WMEQMX— BJ (MMAX)*BY(MMAX- 1) -BJ (MMAX- 1)*BY(MMAX) 
IF(ZABS( (WMEQO  -WEXACT) /WEXACT)  .GT.  1.0D-4  .OR. 


ZABS( ( WMEQMX - WEXACT ) /WEXACT )  .GT.  1.0D-4)then 
write (44,*)'  WARNING -WRONSKIAN  EXCEEDS  TOLERANCE 
write (44 , *) ' ARG  =' ,ARG 

write (44,*) 'WO  =', ZABS( (WMEQO  -WEXACT) /WEXACT) 
write (44 , *) ' WMX= ' , ZABS ( (WMEQMX- WEXACT) /WEXACT) 
end  if 


RETURN 

END 
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SUBROUTINE  CLIN  EQ  (A,  B,  X,  N,  N  DIM,  I FLAG,  ERR) 


C 

C  CLIN  EQ  USES  L-U  DECOMPOSITION  TO  FIND  THE  TRIANGULAR  MATRICES  L 
C  AND  U  SUCH  THAT  L  *  U  -  A.  THE  MATRICES  L  AND  U  ARE  STORED  IN  A. 

C  THIS  FORM  IS  USED  WITH  BACK  SUBSTITUTION  TO  FIND  THE  SOLUTION  X  OF 

C  A*X-L*U*X-B.  NIS  THE  NUMBER  OF  EQUATIONS  AND  NDIM  IS 
C  THE  DIMENSION  OF  ALL  ARRAYS.  ERR  IS  THE  ESTIMATED  RELATIVE  ERROR 
C  OF  THE  SOLUTION  VECTOR. 

C 

C  IF  I FLAG  -  0,  THEN  L,  U  AND  X  ARE  COMPUTED.  OTHERWISE,  IT  IS 

C  ASSUMED  THAT  L  AND  U  HAVE  BEEN  COMPUTED  IN  A  PREVIOUS  CALL  AND  ARE 

C  STILL  STORED  IN  A. 

C 

IMPLICIT  REAL* 8  (A-H,0-Z) 

COMPLEX*16  A,  B,  X,  T 
INTEGER*2  IROW 

DIMENSION  A(N  DIM,  N  DIM),  B(N  DIM),  X(N  DIM) 

DIMENSION  IROW (51) ,  Q(51) 

DATA  EPS  /1.0D-15/ 

C 

IF(N  .GT.  51)  GO  TO  900 
IF (N  .EQ.  1)  GO  TO  850 
IF(IFLAG  .NE.  0)  GO  TO  600 
DO  50  I  -  1 , N 
Q(I)  -  0.0D0 
DO  40  J  -  1,N 
Q Q  -  ZABS(A(I , J) ) 

40  IF(Q(I)  .LT.  QQ)  Q(I)  -  QQ 
IF(Q(I)  .EQ.  0.0D0)  GO  TO  901 
50  CONTINUE 
ERR  -  EPS 
PPIV  -  0.000 
DO  100  I  -  l.N 
100  IROW(I)  -  I 
C 

DO  500  L  —  1 ,N 

PIVOT  -  0.0D0 

K  -  L  -  1 

DO  240  I  -  L,N 

IF(K  .LT.  1)  GO  TO  230 

DO  220  J  -  1 ,K 

220  A(I , L)  -  A(I , L)  -  A(J , L)  *  A(I,J) 


230  F  -  ZABS(A(I ,L) )  /  Q(I) 

IF( PIVOT  .GT.  F)  GO  TO  240 
PIVOT  -  F 
NPIVOT  -  I 
240  CONTINUE 

IF(PIV0T  .EQ.  0.0D0)  GO  TO  901 
IF(PPIV  .LE.  PIVOT)  GO  TO  250 
ERR  -  ERR  *  PPIV  /  PIVOT 
c  PRINT  *, 'ERR-' , ERR 

c  wrlte(44,*) 'err-' ,err 

IF(ERR.GE. 1.0D0)  GO  TO  901 
250  PPIV  -  PIVOT 

IF(NPIVOT  .EQ.  L)  GO  TO  280 
Q(NPIVOT)  -  Q(L) 

J  -  IROW(L) 

IROW(L)  -  IROW(NPIVOT) 


IROW(NPIVOT)  -  J 
DO  260  I  -  1 ,N 
T  -  A(L, I) 

A(L, I)  -  A(NPIVOT , I ) 

A(NPIVOT.I)  -  T 
260  CONTINUE 

280  IF(L  .EQ.  N)  GO  TO  500 
T  -  (1.0D0.0.0D0)  /  A(L,L) 

K  -  L  +  1 

M  -  L  -  1 

DO  450  I  -  K, N 

IF(M  .LT.  1)  GO  TO  400 

DO  350  J  -  1 ,M 

350  A(L, I)  -  A(L, I)  -  A(L, J)  *  A(J,I) 

400  A(L, I)  -  T  *  A(L, I ) 

450  CONTINUE 
500  CONTINUE 

IF(ERR  .GT.  1.0D-5)  PRINT  998,  ERR 

600  DO  620  I  —  2 ,N 
620  X(I)  -  (0 . 0D0 , 0 . 0D0) 

J  -  IROW(l) 

X(l)  -  B(J)  /  A(1 , 1) 

DO  700  I  —  2 ,N 
J  -  IROW(I) 

K  -  I  -  1 
DO  650  L  -  1 ,K 

650  X(I)  -  X(I)  +  A(I , L)  *  X(L) 

X(I)  -  (B( J )  -  X(I) )  /  A(I , I) 

700  CONTINUE 
K  -  N  -  1 
DO  800  I  —  1 ,F 
J  -  N  -  I 
M  -  J  +  1 
DO  800  L  -  M,N 
X(J)  -  X(J)  -  X(L)  *  A(J.L) 

800  CONTINUE 
RETURN 

850  X(1)-=B(1)/A(1,1) 

RETURN 

900  PRINT  999 
ERR  -  1.0D0 
RETURN 

901  PRINT  997 
ERR  -  1.0D0 
RETURN 

997  FORMAT ( ' 0 ERROR  IN  CLIN  EQ ,  MATRIX  IS  SINGULAR') 

998  FORMAT ( ' OCAUTION - ' , 

$  '  CLIN  EQ  HAS  DECOMPOSED  AN  ILL-CONDITIONED  MATRIX.'/ 

$  '  RESULTS  WILL  HAVE  RELATIVE  ERROR  -' .1PE12.2) 

999  FORMAT ( 'OERROR  IN  CLIN  EQ,  MATRIX  SIZE  GREATER  THAN  51') 


It 
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APPENDIX  B 

PLOTTING  ROUTINE  LISTING 


MAIN  ROUTINE  FOR  PROGRAM 


complex*16  cl,c2 

real*8  f rqkhz , deltal , delta2 , xyint , xO , yO , r , alpha 
real*4  xplot .yplotl ,yplot2 ,yplot3 , yplot4 , yplot5 , yplot6 
C0MM0N/C0Ml/xplot(302) ,yplotl(302) ,yplot2(302) ,yplot3(302) , 

yplot4(302) ,yplot5(302) ,yplot6(302) 
COMMON/COM2/frqkhz , deltal , delta2 , cl , c2 , xyint , alpha 
common/com3/iflag,nrslab .numpts 
common/com4/x0 , yO , r 
CHARACTER* 1  ANSWER 
character*20  fn 

call  init 


print  *, 'Enter  name  of  file  containing  data  to  be  plotted: 

read(5,850)  fn 

format(a20) 

OPEN(UNIT-25 , file-fn, status-'OLD' , form-' unformatted' ) 

READ(25)  if lag, nr slab , f rqkhz , deltal , delta2 , cl , c2 , xyint , alpha 
if (if lag  .eq.  0)  then 

READ (25)  xO ,y0 .numpts 
else 

READ(25)  r, numpts 
endif 
icount-0 
do  i— 1, numpts 

read (2 5, end-900)  xplot (i) , yplotl (i) ,yplot2(i) ,yplot3(i) ,yplot4(i) , 
?  yplot5(i) ,yplot6(i) 

xplot (i)-xplot ( i )/l ,0e3 
icount-icount+1 
enddo 


numpts— icount 
CALL  PLTDTA 


SUBROUTINE  PLTBGN 

open(unit-20 , f ile-' /dev/plc7550a' ) 

call  hpinit(2, 0,0, 0,20) 

RETURN 

END 

SUBROUTINE  PLTDTA 
complex*16  cl,c2 

real*8  frqkhz , deltal , delta2 , xyint ,x0 ,y0 , r , alpha 
COMMON/COMl/xplot ( 302 ) , yplotl ( 302 ) ,yplot2(302) ,yplot3(302) , 

$  yplot4(302) ,yplot5(302) ,yplot6(302) 

C0MM0N/C0M2/frqkhz , deltal , delta2 , cl , c2  xyint , alpha 
common/com3/iflag,nr3lab ,numpts 
common/com4/x0 ,y0 , r 

common/com5/  sizex,xmin,xmax,xtic , sizey ,ymin,yniax,ytic 
LOGICAL  UP(302),UPL(2) 

DIMENSION  XL(2) ,YL(2) ,XLAB0(3) ,xlabl(7) , YLABa(5) ,ylabp(6) 
dimension  ia(8) 
character*l  answer 

data(XLAB0(j) ,j-  1,3) /4HRANG , 4HE  (M,4Hm)  / 

data(XLABl(j) , j-  1 , 7)/4HDIST,4HURBA,4HNCE  , 4HL0CA , 4HTION , 4H- S  (, 
$  4HMm)  / 

data(YLABa( j ) , j—  1,5) /4HAMPL, 4HITUD , 4HE  RA , 4HTI0- , 4HdB  / 
data(YLABp(j ) , j-  1,6) /4HPHAS , 4HE  DI , 4HFFER , 4HENCE , 4H  -  D.4HEG  / 

data  ia/82,79,57,48,73,87,73,80/ 

R090IWIP 

HEIGHT-0 . 1 
DO  100  J-1,302 
UP (J)-. FALSE. 

CONTINUE 
DO  110  J-1,2 
UPL(J)-. FALSE. 

CONTINUE 

DRAW  AMPLITUDE  PLOT 

print  *,'Do  you  want  amplitude  plots?' 

read  (5,972)  answer 

format(al) 

if(answer  .eq.  'Y'  .or.  answer  .eq.  'y')  then 
print  * 
print  * 

print  *, 'How  many  amplitude  plots  do  you  want?' 

read  *,  numplt 

print  *, 'numplt-' .numplt 

print  * 

if(numplt  .eq.  1)  then 
print  * 

print  *, 'Enter  appropriate  code  for  desired  plot' 
print  * 

print  *, '  1  ez  plot' 

print  * , '  2  hy  plot' 

print  * , '  3  hr  plot ' 

read  * , iplot 
else 

if (numplt  .eq.  2)  then 


print  * 

print  *, 'Enter  appropriate  code  for  desired  plots' 
print  * 

print  *,'  1  ez  and  hy  plots' 

print  '  2  ez  and  hr  plots' 

print  *,'  3  hy  and  hr  plots' 


print  * , ' 
print  * , ' 
print  * , ' 
read  *,iplot 
end  if 
end  if 


call  input 

CALL  PLTBGN 

call  buff (1 , ia.xbuf , 8) 

CALL  PLOT (1.0, 1.5, -3) 

CALL  BORDER (sizex, xmin, xmax , xtic , 1 , sizey , ymin, ymax,y tic , 1) 

xinc- (xmax -xmin) /sizex 

yinc-(ymax-ymin) /sizey 

XL(1)-1 . 2 

XL(2)-2 . 0 

if((numplt  .eq.  1  .and.  iplot  . eq .  1)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  .eq.  1  .or.  iplot  .eq.  2))  .or. 
$  (numplt  .eq.  3))  then 

CALL  CURVE (xPLOT ,yp loti , UP , numpts ,xmin,ymin,xinc ,yinc , 1) 
CALL  SYMBOL(0. 0,-0. 6, HEIGHT , 8H | EZ/EZU | ,0.0,8) 

YL(1)  — 0.6 
YL(2)--0 . 6 

CALL  CURVE (XL ,YL,UPL, 2, 0.0, 0.0, 1.0,1. 0,1) 
endif 

if((numplt  .eq.  1  .and.  iplot  .eq.  2)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  .eq.  1  .or.  iplot  .eq.  3))  or. 

$  (numplt  .eq.  3))  then 

CALL  CURVE (xPLOT ,yp lot 2 , UP .numpts , xmin , ymin , xinc , yinc ,4) 
CALL  SYMBOL(0.0, -0.8, HEIGHT , 8H | HY/HYU | ,0.0,8) 

YL(1)— -0 . 8 
YL(2)  — 0.8 

CALL  CURVE (XL,YL,UPL,2,0.0,0.0,1.0,1.0,4) 
endif 

if((numplt  .eq.  1  .and.  iplot  .eq.  3)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  . eq .  2  .or.  iplot  .eq.  3))  .or. 

$  (numplt  .eq.  3))  then 

CALL  CURVE (xPLOT , yplot3 ,UP , numpts , xmin , ymin , xinc , yinc , 5) 
CALL  SYMBOL(0 . 0 , - 1 . 0 , HEIGHT , 10H | HR/HRU |  ,0.0,10) 

YL(1)--1.0 
YL(2)-- 1 . 0 

CALL  CURVE (XL, YL.UPL, 2, 0.0, 0.0, 1.0, 1.0, 5) 
endif 

if(iflag  .eq.  0)  then 
xpos-( sizex- 1 . 2)/2 . 0 

CALL  SYMBOL(xpos, -0. 30, HEIGHT, XLAB0.0. ,12) 

else 

xpos-(sizex-2 . 6)/2 . 0 

CALL  SYMBOL(xpos, -0.30, HEIGHT, XLAB1.0. ,28) 
endif 

ypos-( sizey- 1 . 8)/2 . 0 

CALL  SYMB0L(-0. 30, ypos, HEIGHT, YLABa, 90. ,20) 

LABEL  PLOTS 


CALL  SYMBOL( 3.0, -0.6, HEIGHT , 20HFREQUENCY-  HZ ,0.0, 20) 

CALL  NUMBER(4 . 0 , - 0 . 6 , HEIGHT , sngl ( f rqkhz*l . 0e3 ) , 0 . 0 , 1 ) 
if(iflag  .eq.  I)  then 

call  symbol (3.0, -0.8, height, 20hRANGE-  Mm, 0.0, 20) 

call  number (3 .7, -0.8,height,sngl( r/1 . 0e3) ,0.0,3) 
endif 

call  symbol(0 . 0 , sizey+0 . 9 .height , 20hCl=(  ,  ),0.0,20) 

clr-dble(cl) 

cli-dimag(cl) 

CALL  NUMBER (0 . 5 , sizey+0 . 9 .HEIGHT, clr ,0.0,2) 

CALL  NUMBER(1. 3, sizey+0. 9, HEIGHT, cli, 0.0, 2) 

call  symbol (0 . 0 , sizey+0 . 7 .height , 20hC2=(  ,  ),0.0,20) 

c2r-dble(c2) 

c2i-diraag(c2) 

CALL  NUMBER(0. 5, sizey+0. 7, HEIGHT, c2r, 0.0, 2) 

CALL  NUMBER(1. 3, sizey+0. 7, HEIGHT, c2i, 0.0, 2) 
if(iflag  .eq.  0)  then 

call  symbol(0 . 0 , sizey+0 . 5 .height , 14hX0=  KM, 0.0, 14) 

CALL  NUMBER(0. 4, sizey+0. 5, HEIGHT, sngl(xO) ,0.0,2) 

call  symbol ( 3 . 0 , sizey+0 . 5 .height , 14hY0-  KM, 0.0, 14) 

CALL  NUMBER( 3 .4, sizey+0. 5, HEIGHT, sngl (yO) ,0.0,2) 
else 

call  symbol (0 . 0 , sizey+0 . 5 .height , 20hALPHA=  DEG,0.0,20) 

CALL  NUMBER(0 . 7 , sizey+0 . 5 , HEIGHT , sngl (alpha) , 0 . 0 , 2 ) 
if (alpha  .eq.  0.0)  then 

call  symbol(3 . 0 , sizey+0 . 5 .height , 22hY- INTERCEPT*  KM, 

$  0.0,22) 

else 

call  symbol (3 . 0 , sizey+0. 5 .height , 2 2hX- INTERCEPT-  KM, 

$  0.0,22) 

endif 

CALL  NUMBER (4 . 25 , sizey+0 . 5 , HEIGHT , sngl (xyint) ,0.0,2) 
endif 

call  symbol (0 . 0 , sizey+0 . 3 .height , 20hNRSLAB-  ,0.0,20) 

CALL  NUMBER(0. 8, sizey+0. 3, HEIGHT, real (nrslab) ,0.0, -1) 

call  symbol(3 . 0 , sizey+0 . 9 .height , 22hDELTAl«  DEG, 0.0, 22) 

CALL  NUMBER(4. 0, sizey+0. 9, HEIGHT, sngl(deltal) ,0.0,2) 

call  symbol (3 . 0 , sizey+0 , 7 .height , 22hDELTA2-  DEG, 0.0, 22) 

CALL  NUMBER(4 . 0 , sizey+0 . 7 , HEIGHT , sngl (delta2) , 0 . 0 , 2) 

CALL  PLTEND 
endif 


DRAW  SIGNAL  LEVEL  PLOT 

print  *, 'Do  you  want  phase  plots?' 

read  (5,972)  answer 

if(answer  .eq.  'Y'  .or.  answer  .eq.  'y')  then 
print  * 
print  * 

print  *, 'How  many  phase  plots  do  you  want?' 

read  *,  numplt 

print  *, 'numplt-' , numplt 

print  * 

if (numplt  .eq.  1)  then 
print  * 

print  *, 'Enter  appropriate  code  for  desired  plot' 
print  * 

print  '  1  ez  plot' 

print  2  hy  plot' 
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print  *, '  3  hr  plot' 

read  *, iplot 
else 

if (numplt  .eq.  2)  then 
print  * 

print  *, 'Enter  appropriate  code  for  desired  plots' 
print  * 

print  1  ez  and  hy  plots' 

print  *,'  2  ez  and  hr  plots' 

print  3  hy  and  hr  plots' 

read  * , iplot 
end  if 
endif 

call  input 

CALL  PLTBGN 

call  buff (1 , ia.xbuf , 8) 

CALL  PLOT (1.0, 1.5, -3) 

CALL  BORDER(sizex,xmin,xmax,xtic , 1 , sizey,ymin,ymax,ytic , 1) 

xinc-(xmax-xmin) /sizex 

yinc-(ymax-ymin) /sizey 

XL(1)— 1 . 4 

XL(2)-2.2 

if((numplt  .eq.  1  .and.  iplot  .eq.  1)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  .eq.  1  .or.  iplot  .eq.  2))  .or. 
$  (numplt  .eq.  3))  then 

CALL  CURVE (xPL0T,yplot4, UP, numpts ,xmin,ymin, xinc , yinc , 1) 
CALL  SYMBOL(0 . 0 , -0 . 6 , HEIGHT , 13HPHASE(EZU/EZ) ,0.0,13) 

YL(1)  —  0.6 
YL(2)  — 0.6 

CALL  CURVE (XL , YL , UPL ,2, 0.0, 0.0, 1.0, 1.0,1) 
endif 

if ((numplt  .eq.  1  .and.  iplot  . eq .  2)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  .eq.  1  .or.  iplot  .eq.  3))  .or. 

$  (numplt  .eq.  3))  then 

CALL  CURVE(xPL0T,yplot5 , UP .numpts ,xmin ,ymin , xinc , yinc ,4) 
CALL  SYMBOL(0 . 0 , -0 . 8 , HEIGHT , 13HPHASE(HYU/HY) ,0.0,13) 

YL(1)  —  0.8 
YL(2)~  0.8 

CALL  CURVE (XL, YL, UPL, 2,0. 0,0. 0,1.0, 1.0,4) 
endif 

if ((numplt  .eq.  1  .and.  iplot  .eq.  3)  .or. 

$  (numplt  .eq.  2  .and.  (iplot  .eq.  2  .or.  iplot  .eq.  3))  .or. 

$  (numplt  .eq.  3))  then 

CALL  CURVE (xPLOT ,yplot6 , UP .numpts , xmin,ymin, xinc , yinc , 5) 
CALL  SYMBOL(0.0, -1.0, HEIGHT, 13HPHASE(HRU/HR) ,0.0,13) 

YL(1)  — 1.0 
YL(2)  — 1.0 

CALL  CURVE ( XL, YL, UPL, 2, 0.0, 0.0, 1.0, 1.0, 5) 
endif 

if (if lag  .eq.  0)  then 
xpos-(sizex-l . 2)/2 . 0 

CALL  SYMBOL (xpos ,-0.30, HEIGHT , XLAB0 , 0 . , 12 ) 
else 

xpos-(sizex-2 . 6)/2.0 

CALL  SYMBOL(xpos, -0.30, HEIGHT, XLAB1.0. ,28) 
endif 


ypos-(sizey-2 . l)/2 . 0 

CALL  SYMBOL( - 0 . 30 , ypos , HEIGHT , YLABp , 90 . , 24 ) 

LABEL  PLOTS 

CALL  SYMBOL( 3 . 0 , -0 . 6 , HEIGHT , 20HFREQUENCY-  HZ , 0 . 0 , 20) 

CALL  NUMBER (4 . 0 , - 0 . 6 , HEIGHT , sngl ( f rqkhz*l . 0e3 ) , 0 . 0 , 1 ) 
if(iflag  .eq.  1)  then 

call  symbol (3 . 0 , -0 . 8 .height , 20hRANGE—  Mm, 0.0, 20) 

call  number (3 . 7 , -0. 8 .height , sngl(r/l . 0e3) , 0. 0, 3) 
endif 

call  symbol (0 . 0 , sizey+0 . 9 .height , 20hCl-(  ,  ),0.0,20) 

clr-dble(cl) 

cli-dimag(cl) 

CALL  NUMBER(0 . 5 , sizey+0 . 9 , HEIGHT , clr , 0 . 0 , 2) 

CALL  NUMBER (1 . 3 , sizey+0 . 9 , HEIGHT , cl i , 0 . 0 , 2 ) 

call  symbol(0 . 0 , sizey+0 . 7 .height , 20hC2-(  ,  ),0.0,20) 

c2r-dble(c2) 

c2i-dimag(c2) 

CALL  NUMBER (0 . 5 , sizey+0 . 7 , HEIGHT, c2r , 0.0,2) 

CALL  NUMBER ( 1 . 3 , sizey+0 . 7 , HEIGHT , c2 i , 0 . 0 , 2 ) 
if (if lag  .eq.  0)  then 

call  symbol(0 . 0 , sizey+0 . 5 .height , 14hX0“  KM, 0.0, 14) 

CALL  NUMBER(0. 4, sizey+0. 5, HEIGHT, sngl (xO) ,0.0,2) 

call  symbol (3 . 0, sizey+0. 5 .height , 14hY0=  KM, 0.0, 14) 

CALL  NUMBER( 3. 4, sizey+0. 5, HEIGHT, sngl (yO) ,0.0,2) 
else 

call  symbol (0 . 0 , sizey+0 . 5 .height , 20hALPHA=  DEG, 0.0, 20) 

CALL  NUMBER ( 0 . 7 , s izey+0 . 5 , HEIGHT , sngl ( alpha ) , 0 . 0 , 2 ) 
if(alpha  .eq.  0.0)  then 

call  symbol(3.0, sizey+0. 5, height, 2 2hY- INTERCEPT-  KM, 

$  0.0,22) 

6 1 S  0 

call  symbol (3 . 0 , sizey+0 . 5 , height , 22hX- INTERCEPT*  KM, 

$  0.0,22) 

endif 

CALL  NUMBER (4 . 25 , sizey+0 . 5 , HEIGHT , sngl (xyint) ,0.0,2) 
endif 

call  symbol (0 . 0 ,  sizey+0 . 3  .height , 20hNRSLAB*=  ,0.0,20) 

CALL  NUMBER(0 . 8 , sizey+0 . 3 , HEIGHT, real (nrslab) , 0 . 0 , - 1) 
call  symbol(3 . 0 , sizey+0 . 9 .height , 22hDELTAl-  DEG, 0.0, 22) 

CALL  NUMBER (4 . 0 , sizey+0 . 9 , HEIGHT, sngl (deltal) ,0.0,2) 
call  symbol ( 3 . 0 , sizey+0 . 7 .height , 22hDELTA2=  DEG, 0.0, 22) 

CALL  NUMBER (4 . 0 , sizey+0 . 7 , HEIGHT , sngl(delta2) ,0.0,2) 
r''  "LTEND 
;if 

KETUR 

END 

SUBROUTINE  P LTEND 
CALL  NEWPEN(O) 

CALL  PLOT ( 0 . ,0. ,999) 

RETURN 

END 


subroutine  input 

common/com5/  sizex ,xmin,xmax ,xtic , sizey , ymin, ymax, y tic 
character*l  answer 

print  * 

print  *, 'Current  values  for  sizex  and  sizey  are sizex, sizey 
print  *,'Do  you  want  to  change  them?' 
read  (5,972)  answer 
972  format(al) 

if(answer  .eq.  'y'  .or.  answer  .eq.  'Y')  then 
print  *, 'Enter  values  for  sizex, sizey' 
read  *, sizex, sizey 
print  *,' size-' , sizex, sizey 
endif 

print  * 

print  *, 'Current  values  for  xmin  and  xmax  are: ' ,xmin,xmax 
print  *,'Do  you  want  to  change  them?' 
read  (5,972)  answer 

if (answer  .eq.  'y'  .or.  answer  .eq.  'Y')  then 
print  *, 'Enter  values  for  xmin, xmax' 
read  *, xmin, xmax 
print  *, 'xmin, xmax-' , xmin, xmax 
endif 

print  * 

print  *, 'Current  values  for  ymin  and  ymax  are : ' ,ymin,ymax 
print  *,'Do  you  want  to  change  them?' 
read  (5,972)  answer 

if (answer  .eq.  'y'  .or.  enswer  .eq.  'Y')  then 
print  *, 'Enter  values  for  ymin, ymax' 
read  *, ymin, ymax 
print  *,' ymin , ymax- ', ymin , ymax 
endif 

print  * 

print  *,'Tic  marks  will  be  every  ',xtic,'  units  on  the  x-axis  and' 
print  *, '  every  ',ytic,'  units  on  the  y-axis.' 

print  *,'Do  you  want  to  change  them?’ 
read  (5,972)  answer 

if (answer  .eq.  'y'  .or.  answer  .eq.  'Y')  then 
print  *, 'Enter  values  for  xtic.ytic' 
read  *, xtic.ytic 
print  *, 'xtic.ytic-' , xtic.ytic 
endif 
return 
end 


c 


subroutine  border (xlng.xmin, xmax, xinc , nx, ylng, ymin, ymax, yinc ,ny) 


dimension  xinc(nx) ,yinc(ny) 
logical  fy.fx 

fx= . false . 
fy= . false . 

if(nx  .eq.  1)  fx=.true. 
if(ny  .eq.  1)  fy-.true. 
xt=xlng- . 1 
yt-ylng- .1 

xscale=xlng/(xmax-xmin) 
yscale=ylng/( ymax -ymin) 
ym-abs (ymin) 
yin-- . 4 

if(ym  .ge.  10.)  y ln-yln-. 1 
iffym  . gt .  ICO.)  yin— yin-. 1 
if(ym  . ge .  1000.)  yln=yln- . 1 
if(ymin  .It.  0.)  yln=yln- . 1 
ym— abs (ymax) 
ylm-- . 4 

if(ym  . ge .  10.)  ylm=ylm-.l 
if(ym  .ge.  100.)  ylm=ylm- . 1 
if(ym  .ge.  1000.)  ylm-ylm-.l 
if (ymax  .It.  0.)  ylm-ylm-.l 
xm-abs (xmax) 
xlm=- . 3 

if(xm  . ge .  10.)  xlm-xlm- . 1 
if(xm  . ge .  100.)  xlm=xlm-.l 
if (xm  .ge.  1000.)  xlm=xlm-.l 
if (xmax  .It.  0.)  xlm-xlm-.l 
if(fx)  dx-xinc(l) 
if(fy)  dy-yinc(l) 

iy-i 

yi-°. 

call  number (yin , 0 .,. 1 ,ymin, 0 ., 1) 
call  plot (0 . , 0 . ,3) 
if(fy)  go  to  110 

10  yp-(yinc(iy) -ymin)*yscale 
go  to  111 

110  yl-yl+dy 
yp-yl*yscale 

111  if (yp  .It.  0.)  go  to  99 
if (yp  .ge.  ylng)  go  to  11 
call  plot(0 . ,yp , 2) 

call  plot( . 1 ,yp , 2) 
call  plot(0 . ,yp , 2) 
if(fy)  go  to  110 
iy-iy+l 

if(iy  . le .  ny)  go  to  10 

11  call  plot(0 . ,ylng, 2) 

call  number (ylm, ylng- . 1 , . 1 ,ymax ,0 . , 1) 
call  plot(0. ,ylng, 3) 
ix-1 
xl-0 . 

if(fx)  go  to  112 

12  xp-(xinc(ix) -xmin)*xscale 
go  to  120 


120 


112  xl-xl+dx 
xp-xl*xscale 

120  if(xp  .It.  0.)  go  to  99 
if(xp  .ge.  xlng)  go  to  13 
call  plot(xp ,ylng, 2) 
call  plot(xp ,yt , 2) 
call  plot(xp , ylng, 2) 
lf(fx)  go  to  112 
ix-ix+1 

if(ix  . le .  nx)  go  to  12 

13  call  plot(xlng,ylng, 2) 
if(fy)  go  to  130 

113  iy-iy-1 

if(iy  .le.  0)  go  to  15 
yp-(yinc(iy) -ymin)*yscale 
go  to  14 

130  yl-yl-dy 

yp-yl*y scalp 

if(yp  .le.  0.)  go  to  15 

14  call  plot(xlng,yp , 2) 
call  plot(xt ,yp , 2) 
call  plot(xlng,yp, 2) 
if(fy)  go  to  130 

go  to  113 

15  call  plot(xlng,0. ,2) 

call  number (xlng+xlm, - . 2 , . 1 ,xmax ,0 . , 1) 
call  plot(xlng, 0 . , 3) 
if(fx)  go  to  150 

115  ix-ix-1 

if(ix  .le.  0)  go  to  17 
xp-(xinc(ix) -xmin)*xscale 
go  to  16 

150  xl-xl-dx 

xp-xl*xscale 

if(xp  .le.  0.)  go  to  17 

16  call  plot(xp , 0 . , 2) 
call  plot(xp , . 1 , 2) 
call  plot(xp,0. ,2) 
if(fx)  go  to  150 
go  to  115 

17  call  plot(0. ,0. ,2) 

call  number (0 2 ,. 1 ,xmin, 0 1) 
return 

99  print  100 , xlng, xmin, xmax, xinc ( 1) ,nx, ylng, ymin,ymax,yinc( 1) ,ny 

100  format (’0***  Error  in  BORDER:  xlng,  xmin,  xmax,  xinc(l),  nx  , 

$  Ip4el5 . 5 , i5/24x, 'ylng ,  ymin,  ymax,  yinc(l),  ny  -’,lp4el5.5, 

$  i5/'0***') 

call  pltend 

stop 

end 
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subroutine  curve (x,y , up , nrpts , xmin,ymin,xinc ,yinc , line) 


* 
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x.y.up  must  be  dimensioned  at  least  nrpts 
xmin.ymin  are  x,y  origin  in  user  units 
xinc.yinc  are  x,y  scales  in  user  units  per  inch 


line-1 

2 


solid 
long  dash 
medium  dash 
short  dash 
dotted 

short  +  long  dash 
short  +  short  +  long  dash 


logical  up,upl,up2 

dimension  ipen(8) ,  joc(7) ,x(nrpts) ,y(nrpts) ,up(nrpts) 
data  ipen/2 , 2 , 2 , 3 , 2 , 3 , 2 , 3/ , j oc/18 ,  11,  14,  23,  32,  41,  16/ 
data  delr/.l/ 


if (nrpts  .le.  1)  go  to  99 


if (line)  1,2,3 

kk-mod(line , 7) +7 

go  to  4 

kk-0 

go  to  4 

kk=mod(line , 7) 

kk=kk+l 

jo=joc(kk)/10 

j  c=j  oc (kk) - 10*j  o 

ip-ipen(jo) 


j“0 
dr=0 . 
rhol=0 . 
rho2=delr 

pxl=(x(l) -xmin)/xinc 
pyl-(y(l) -ymin)/yinc 
upl— up ( 1 ) 

if (.not.  upl)  then 


go  to  first  position  with  pen  up 
call  plot(pxl ,pyl , 3) 
if(kk  ,eq.  6)  then 
px2-(x(2) -xmin)/xinc 
py2-(y(2) -ymin)/yinc 
delx-px2-pxl 
dely-py2-pyl 

rho-sqr t ( delx**2+dely**2 ) 
if(rho  .eq.  0.)  then 
dx  6-delx*. 1 
dy  6-dely*.l 
else 

dx  6-delx*delr/rho*. 1 
dy  6-dely*delr/rho* . 1 
end  if 

call  plot(pxl+dx6 ,pyl+dy6 , 2) 
end  if 


do  40  i-2,nrpts 
px2-.(x(i)  -xmin)/xinc 
py2-(y(i)  -ymin)/yinc 
up2-up(i) 
i£(up2)  then 
dr-0. 
rhol-0. 
rho2-delr 
go  to  39 
end  if 

if(upl)  then 

c  pen  has  been  up,  prepare  to  lower  pen 
call  plot(px2 ,py2 , 3) 
go  to  39 
end  if 

if(kk  .eq.  2)  go  to  38 

delx-px2-pxl 

dely-py2-pyl 

rho-sqrt(delx**2+dely**2) 

rhol-rhol+rho 

if(rho2  .gt.  rhol)  go  to  38 

delx-delx*delr/rho 

de ly-de ly*de 1 r/rho 

dx  6-delx*.l 

dy  6— dely*.l 

if (dr  .eq.  0.)  go  to  20 

dx-de lx*dr /de 1 r 

dy-dely*dr/delr 

pxl-pxl+dx 

pyl-pyl+dy 

go  to  21 

20  if(rho2  .gt.  rhol)  go  to  38 

pxl-pxl+delx 
pyl-pyl+dely 

21  call  plot(pxl ,pyl , ip) 

if(kk  .eq.  6)  call  plot(pxl+dx6 ,pyl+dy6 , 2) 
j-j+l 

ip-ipen( jo+mod( j , jc)) 
rho2-rho2+delr 
go  to  20 

38  call  plot(px2 ,py2 , ip) 
dr-rho2-rhol 

39  pxl-px2 
pyl-pyt 
upl-up2 

40  continue 

99  return 

end 


